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Chapter 1 


Introduction 


Both the practitioner using computational fluid dynamics in engineering design 
calculations and the scientist grappling with gaps in the theoretical foundations 
are aware that much remains to be done before the subject can be put on 
firm ground. This is particularly true in the theory and numerical analysis 
of hyperbolic conservation laws, vital in gas dynamics and compressible fluid 
mechanics and a fundamental component in the resolution of solutions of the 
Navier-Stokes equations for compressible flow. There the numerical solution of 
hyperbolic systems is confronted with a list of major difficulties and questions 
that have been under study for many years. 

These include classical problems of numerically resolving shocks and dis- 
continuities, characteristic of solutions of hyperbolic problems, while simulta- 
neously producing high-order, non-oscillatory results near shocks and elsewhere 
in the solution domain. Moreover, the basic issue of quality of numerical solu- 
tions is fundamentally important: how accurate are the numerical simulations 
and how does one obtain the most accurate results for a fixed computational re- 
source? These questions lie at the core of modern adaptive methods that aim 
to control the error in the computed solution and to optimize the computa- 
tional process. In addition, methodologies that attempt to address these issues 


1 
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cannot be limited to one-dimensional cases; they must be extendable to prob- 
lems involving realistic geometries, boundary and initial conditions in arbitrary 
domains in two- and three-dimensions. Finally, there is the issue of computa- 
tional efficiency. Modern numerical schemes for large-scale applications should 
be readily parallelizable for implementation in emerging multi-processor archi- 
tectures. 

This dissertation addresses these issues for model classes of hyperbolic 
conservation laws. The basic approach developed in this work employs a new 
family of adaptive, /ip-version, finite element methods based on a special dis- 
continuous Galerkin formulation for hyperbolic problems. The discontinuous 
Galerkin formulation admits high-order local approximations on domains of 
quite general geometry, while providing a natural framework for finite element 
approximations and for theoretical developments. The use of /ip-versions of the 
finite element method makes possible exponentially convergent schemes with 
very high accuracies in certain cases; the use of adaptive hp - schemes allows 
/i-refinement in regions of low regularity and p-enrichment to deliver high ac- 
curacy, while keeping problem sizes manageable and dramatically smaller than 
many conventional approaches. The use of discontinuous Galerkin methods is 
uncommon in applications, but the methods rest on a reasonable mathemat- 
ical basis for low-order cases and has local approximation features that can 
be exploited to produce very efficient schemes, especially in a parallel, multi- 
processor environment. 

The place of this work is to first and primarily focus on a model class 
of linear hyperbolic conservation laws for which concrete mathematical re- 
sults, methodologies, error estimates, convergence criteria, and parallel adap- 
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tive strategies can be developed, and to then briefly explore some extensions to 
more general cases. Next, we provide preliminaries to the study and a review 
of some aspects of the theory of hyperbolic conservation laws. We also provide 
a review of relevant literature on this subject and on the numerical analysis of 
these types of problems. 


1.1 Some Mathematical Preliminaries 


The broad aim of this work is to lay mathematical foundations and to develop 
new numerical schemes which will aid in understanding and solving general 
systems of hyperbolic conservation laws of the form 

U >t + F‘(U) iari = S,x = {x l ,x 2 ) e ft C $ d ,t > 0 (1.1) 

where U is a vector consisting of m components to be conserved in fl and 
S is a vector-valued source term. The flux vectors F *,* = 1 , •••,</ are, in 
general, a nonlinear function of U. The subscript t denotes differentiation with 
respect to time and the subscripts X{ denote differentiation with respect to the 
spatial coordinate a?,-, and in which summation convention for repeated indices 
is employed. Alternatively, (1.1) can be written in the quasi-linear form 

U,t + A’(U)U >a;< = S (1.2) 


where the flux Jacobian matrices 


A ! = 


dF i 

dU 


(1.3) 


have real eigenvalues. To complete the initial-boundary value problem, one 
must specify initial conditions of the form 


U(x,0) = U o (x) 


(1.4) 
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and appropriate boundary conditions. 

Most of this investigation focuses on a simpler and restricted class of 
hyperbolic problems for which concrete results can be obtained and for which it 
is possible to treat some specific numerical issues in more detail than is possible 
for the general case. 

The first difficulty encountered in developing numerical methods for 
hyperbolic systems of conservation laws in multiple space dimensions is the 
absence of a general mathematical theory. Hence, many popular numerical 
schemes for multi-dimensional systems are based on existence and uniqueness 
results for scalar hyperbolic conservation laws of the form 

u,t + r(u),xi = 0; x € t > 0 (1.5) 

u(x,0) = u 0 (x) 

where the fluxes /', i = 1, • • • , d are Lipschitz continuous functions of u and 
the initial data no £ L°°(^R d ) fl Z/ 1 (3? d ) have compact support. 

Historically, existence and uniqueness proofs have relied upon compact- 
ness arguments for sequences of solutions generated by the vanishing viscos- 
ity method, see Kruzkov [32], or low-order finite difference approximations, 
see Glimm [22] and Crandall and Majda [15]. More recently, uniqueness re- 
sults have been generalized using the concept of measure- valued solutions (see 
DiPerna [20]) providing a new tool for convergence proofs for a variety of numer- 
ical methods [14], [31]. Here some well-known results for scalar conservation 
laws are summarized. 

Solutions to (1.5) develop discontinuities in finite time, even if the initial 
data is smooth. Thus, (1.5) must be interpreted in the sense of distributions, 



or equivalently, weak solutions are sought which satisfy 

/ / (u<j> t + ^2 f'(u)<f>xi] dx. dt + f u 0 (j> dx ~ 0 (1.6) 

JO J \ ) Jftd 

for all test functions <f> £ C°°(3? d x [0, T )) with compact support. There are an 
infinite number of weak solutions satisfying (1.6) for given initial data u 0 ; how- 
ever, the physically relevant solution satisfies an additional constraint, namely 
an entropy condition. The entropy condition takes into account the fact that 
physical processes are dissipative and that (1.5) models a physical process in 
the limit as the dissipation tends to zero. This solution, the so-called entropy 
solution, satisfies 

f f \u - c\<f) t + sign(u - c)(f l (u) - f(c))(j) Xi dx dt >Q (1.7) 
J 0 J 3 ?^ 

for all nonnegative test functions f x [0,T)) and all c £ :R. 

The following lemma summarizes some basic properties of solutions to 
(1.5) or (1.6). 

Lemma 1 (Crandall and Majda [15]) For every choice of initial data uq £ 
L°°(3ff rf )n L 1 (3f? <i ) ; there exists a unique entropy solution u £ (7([0,oo) : L 
of (1.5) with tt(x,0) = uo(x). Denoting u(x, t) by E(i)v,Q, we have 

(i) \\E(t)uo — ^(0 U o||n(»<*) < ||u 0 — 1 1 £,1 ( 3 J<i) 

(ii) u 0 < v 0 a.e. — » E(t)u 0 < E(t)v 0 a.e. 

(in) uq £ [a, 6] a.e. — > E(t)uo £ [a, b] a.e. 

(iv) Ifuo £ BV(?R d ), t — > E(t)uo is Lipschitz continuous into T 1 (3? li ) 
and ||2?(*) u o||bv(R' 1 ) < H^ollBvq^d) 
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The functions in the space BV($t?) C Lj oc (^R d ) have bounded variation 
and distributional derivatives that are locally measures. The variation of a 
function v 6 BV(dt d ) for d = 2 is defined as 

«ese\'{o}^3i 2 |£| 

|u(x,j/ + ^) -v(x,y ) | 


kllev(^) = sup J 


+ sup 

5 €»\{ 0}-^ 2 


J& 




-dx dy 
\-dx dy 


These results are typical of the mathematical theory that has influenced 
the development of many numerical schemes over the past decade: they apply 
to unbounded domains, scalar-valued solutions, and characterize the solutions 
as among function classes of low regularity. Our goal is to consider problems 
in bounded domains with specific inflow conditions, since these are the types 
of problems encountered in realistic application of the theory. Moreover, when 
there exist portions of the domain over which the solution is smooth, we wish to 
take advantage of that smoothness by exploiting higher-order methods which 
can exhibit high local accuracies. These types of considerations suggest discon 
tinuous Galerkin methods as an approach worthy of study and set the present 
work apart from most conventional approaches to this subject. 


1.2 Higher-Order Methods 

A major challenge in designing higher-order methods for the numerical solu- 
tion of hyperbolic conservation laws over the last decade has been to prevent 
nonphysical oscillations from occuring near discontinuities without destroying 
accuracy in smooth regions. These oscillations can pollute the solution glob- 
ally and can lead to numerical instabilities not revealed by traditional stability 
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analysis applied to a linearized equation. One should also note that the notion 
of the ’’order” of a numerical scheme is not always defined consistently in the 
literature. In classical finite-difference literature, the order of a scheme refers 
to the order of the truncation error in the time-step At, and this may be quite 
different from the actual order of the error in, say, the L 2 -norm. In discussing 
the order of various schemes in the literature on this subject, we generally refer 
to the truncation error, since this use is prevelent, albeit imprecise. In subse- 
quent Chapters, we develop error estimates in well-defined norms so that the 
question of order of accuracy is clearly resolved. 

The classical remedy for controlling oscillations, namely, regularizing 
the conservation law by adding an ’’artificial” diffusive term (e.g. [28], [34], 
and [34]), is easily applied to methods of arbitrarily high order [18], [49]. Un- 
fortunately, this approach is not completely effective at eliminating oscillations 
and may destroy accuracy in smooth regions. The most successful methods 
for solving realistic problems are typically second-order accurate and attack 
oscillations directly by simply preventing them from occuring. These schemes, 
such as those based on the Flux- Corrected Transport (FCT) ideas of Boris 
and Book [8], or the Total Variation Diminishing (TVD) ideas of Harten [25], 
and the monotone reconstruction ideas of Van Leer [51], use a form of flux 
or solution limiting to enforce the discrete counterpart of properties (iii) or 
(iv) in Lemma 1 onto the approximate solution. Unfortunately, this limiting is 
based on one-dimensional concepts which are heuristically extended to multi- 
dimensional systems. Often these extensions result in a loss of accuracy in 
smooth regions [24]. 

Higher-order accurate (greater than second-order) methods for discon- 
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tinuous solutions of hyperbolic conservation laws are primarily in a develop- 
mental stage. Spectral methods have been combined with FCT ideas [10] and 
filtering methods [50] to control oscillations, but much of the work is for one 
space dimension and lacks the geometric flexibilty needed for adaptivity. Es- 
sentially NonOscillatory (ENO) schemes introduced by Harten, Enquist, Os- 
her and Chakravarthy [26] produce high-order approximations for hyperbolic 
problems by using high-order polynomial interpolation of solution mean values. 
Oscillations are controlled by using a solution dependent stencil which avoids 
interpolation across discontinuities. To date, a full theoretical basis for these 
schemes is not available beyond one dimensional cases. 

Cockburn, Shu, and collaborators [12], [13] developed one of the first 
high-order numerical schemes for hyperbolic conservation laws in two space 
dimensions. This work employed Runge-Kutta schemes for advancing the so- 
lution in time and an elaborate projection strategy that guaranteed that the 
total variation of the solutions remained bounded throughout the evolution pro- 
cess. Their TVB (total variation bounded) schemes thus generalized the TVD 
schemes of Harten and others. Goodman and LeVeque [24] showed that TVD 
schemes are at most first-order accurate in dimensions greater than one and 
hence, are not justified mathematically in two- or more dimensions. The TVB 
schemes, however, provide a basis for the development of high-order schemes 
on spatial domains of dimensions two and three. The Cockburn and Shu TVB 
methods are constructed using discontinuous Galerkin methods and, thus, ex- 
tend these methods to nonlinear conservation laws. However, the emphasis of 
their work was in treating the discontinuous Galerkin method as a higher-order 
finite volume method, that is, focusing on the accuracy of solution mean values. 



Moreover, these TVB schemes, thusfar, have delivered only second-order meth- 
ods on non-cartesian meshes in multiple space dimensions. This work inspired 
some of the developments discussed in Chapter 7. 

Among the earliest work on finite element approximations of hyperbolic 
problems is the classical paper of Lesaint and Raviart [33] which introduced 
the discontinuous Galerkin method for linear scalar hyperbolic problems. This 
work contained the first a priori error estimates for h-ve rsion methods based 
on elements of arbitrary, but uniform, polynomial order p . In their work, sub- 
optimal error estimates, with a loss in global accuracy of 0(h) in the L 2 -norm, 
were obtained. 

A detailed analysis of discontinuous Galerkin methods for h-ve rsion 
methods was contributed by Johnson and his collaborators [29], [30]. There 
quasi-optimal a priori estimates showed a global accuracy of 0(/i p+ 2 ) in mesh- 
dependent norms. This work provided a general approach to the mathematical 
analysis of these methods that proved to be invaluable in the present work. 
Among the results established in the present study are developments of new 
a priori and a posteriori error estimates for /ip-version discontinuous Galerkin 
finite element approximations of linear, scalar hyperbolic conservation laws. 
Thus, this study extends and generalizes the results of Johnson and others to 
p- and hp~ve rsion finite elements and provides, for the first time, a posteri- 
ori error estimates for such problems using extensions of the element residual 
method (see [1], [2]). 

One reason for renewed interest in discontinuous Galerkin methods is 
the advent of parallel computations. The decomposition of large-scale prob- 
lems into several computational components that can be handled simultane- 
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ously by multiple processors makes possible significant improvements in the 
efficiency with which large hyperbolic systems can be resolved. Some progress 
in parallelizations of high-order schemes for hyperbolic problems in one- and 
two-dimensions has been made by Biswas, Devine, and Flaherty [7]. in their 
work, extensions of the ideas of Cockburn and Shu [12], [13] are presented 
which make use of higher moments of the solutions over an element in defining 
a broader class of projections for imposing TVB behavior on the entire solu- 
tion in an element and not just the mean values. At present, however, moment 
limiting is restricted to cartesian grids. A component of the work reported 
here is concerned with parallel computing methods for ftp-finite element ap- 
proximations of scalar conservation laws. The fact that discontinuous Galerkin 
methods involve very localized approximations over individual elements makes 
our techniques particularly amenable to element -by-element decomposition and 
parallel processing. 

To a great extent, the present work represents a significant departure 
from conventional methods for the numerical solution of hyperbolic problems. 
Several fundamental issues are addressed: the use of discontinuous ftp-methods, 
to provide high-order local approximation to deliver high accuracy when possi- 
ble but also allowing mesh refinement to resolve irregularities in the solutions; 
the development of a priori error estimates to establish proofs of convergence 
and qualitative information on the performance of the method; the development 
of a posteriori error estimates to monitor the performance of the calculation and 
to estimate quality of the solution; the development of new adaptive strategies 
to control error and optimize meshes; and the development of parallel comput- 
ing strategies to exploit the local character of discontinuous approximations 



and to increase the speed with which solutions can be obtained. In addition 
to these developments, applications to nonlinear conservation laws are also 
presented. 

1.3 A Posteriori Error Estimation and Adaptivity 

The power of adaptivity to efficiently improve solution accuracy was recognized 
early on in the development of unstruct ured grid methods for hyperbolic conser- 
vation laws. These ^-adaptive methods, based on refinement and derefinement 
of an initial mesh [35], [36], [16] or a complete remeshing of the domain [43], 
continue to be the preference for realistic flow simulations. With an emphasis 
on resolving certain features of the solution, many refinement indicators have 
been proposed which are based on some a priori knowledge of the solution 
behavior associated with certain phenomena. Typically, these indicators are 
loosely based on interpolation error estimates applied to key variables and can 
be grossly in error. While this approach may provide some relative measure of 
the local error in the solution, it does not in general provide a reliable estimate 
of the actual error in the approximate solution. 

While the bulk of previous work has concentrated on ^-adaptive meth- 
ods combined with low-order approximations, the effectiveness of p-adaptive 
[7] and ftp-adaptive [18] methods has been demonstrated for certain classes 
of hyperbolic problems. As noted earlier, the present study extends existing 
adaptive strategies to discontinuous approximations of hyperbolic conservation 
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1.4 Scope 

Following this introduction, a new formulation of the discontinuous Galerkin 
method is given for a model class of steady-state, scalar, linear hyperbolic prob- 
lems in two dimensions. There a notion of ftp-dependent norms is introduced 
which generalizes to ftp-methods the idea of mesh-dependent norms used by 
Johnson and Pitkaranta [30], Conceptually, one considers a partition of a do- 
main C J? 2 into finite elements and assigns to each element K a positive 
number p K which is designed to appear in coefficients of a mesh-dependent 
norm in a way to optimize subsequent estimated convergence rates. The num- 
bers p K are identified with the maximum spectral orders of the shape functions 
used in approximations over K, A priori error estimates are derived in these 
norms. 

In Chapter 3, the subject of a posteriori error estimates for the model 
problem is investigated. An extension of the element residual method to hy- 
perbolic conservation laws is described. In the present investigation, two types 
of estimates are produced, one which delivers an upper bound to the global 
error in a suitable mesh-dependent norm and a lower bound in another related 
norm. Theorems are proven which establish that these estimates are indeed 
valid bounds on appropriate measures of the approximation error. 

The availability of both a priori and a posteriori error estimates provides 
a powerful basis for developing adaptive strategies to control the error. In 
Chapter 4, an extension of the 3-step adaptive strategy for ftp-finite element 
methods is presented. This work extends the development in [39] to hyperbolic 
conservation laws and represents, we believe, the first ftp-adaptive methodology 
ever developed for this class of problem. 
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Chapter 5 is devoted to numerical experiments and testing of the the- 
oretical results developed in earlier chapters. Several model problems in two 
dimensions are studied, including examples in non-rectangular domains. The 
numerical results exhibit significant features of the theory and the methodolo- 
gies developed: 1) the asymptotic rates of convergence predicted by our theory 
of a priori estimates are fully confirmed by the computed rates; 2) exponential 
rates of convergence or super algebraic rates are observed, justifying finally the 
decision to use nonuniform hp - meshes for these types of problems; 3) the a 
posteriori estimation methods produce very good estimates of the actual error, 
with effectivity indices near unity in many cases, and with remarkably good 
local error indicators in most of the cases considered; 4) the 3-step adaptive 
strategy works surprisingly well and delivers a targeted error level quite regu- 
larly in around 3 steps; the significant accuracy with which the 3-step scheme 
was able to produce a mesh with a prescribed global error was unexpected. 

In Chapter 6, issues of parallelization of the adaptive methods are in- 
vestigated. There the local character of discontinuous Galerkin methods is 
exploited when possible. A parallel algorithm designed for implementation on 
the Intel iPSC860 computer with 16 processors is described and the results of 
numerical tests are presented. A major issue in the parallel implementation 
of /ip-adaptive schemes is the design of domain decomposition strategies which 
result in a balanced work load on all processors. Two domain decomposition 
strategies recently developed by Patra [40] are used and nearly linear increases 
in speed are observed for certain cases. 

Chapter 7 is devoted to some preliminary results on extensions of the 
work to nonlinear conservation laws. A model problem involving the solu- 
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tion of Burgers’ equation on a two-dimensional domain is investigated. The 
exploratory results suggest that the methods developed here could be very ef- 
fective for problems of this type. The local projection developed in this work 
is extremely effective at controlling oscillations at discontinuities without de- 
stroying accuracy in smooth regions. The projection strategy is a simple one 
designed with the idea of combining low-order approximations at shocks with 
higher-order approximations in smooth regions. 

Finally, in Chapter 8, the major conclusions of the study are given 
together with suggestions for future work. 



Chapter 2 


The Discontinuous Galerkin Method 


The methods presented in this chapter are valid for hyperbolic systems of con- 
servation laws in multiple space dimensions. For clarity of presentation and 
for the purposes of analysis, however, we limit the discussion to a scalar lin- 
ear hyperbolic conservation law. We begin with a detailed description of the 
method for a linear model problem and prove some important properties. Next 
we describe a finite element approximation and derive a priori error estimates 
for an hp-ve rsion of the discontinuous Galerkin method. 

2.1 A Linear Model Problem 

We consider a linear scalar hyperbolic conservation law on a convex polygonal 
domain fl. Let = (/?i,/ 3 2 ) T denote a constant unit velocity vector. The 
domain boundary dQ, with an outward unit normal vector n(x) consists of two 
parts: an inflow boundary T_ = {x € <9fl | /3 • n(x) < 0} and an outflow 
boundary T + = dfl \ T_. Let u denote the quantity that is to be conserved in 
fl and consider the following hyperbolic boundary-value problem: 


/3 • Vii + au = f in 11 C 71 2 


(2.1) 


15 
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(3 ■ n u = (3 ■ n g on T_ (2-2) 

where / G L 2 (Q), g G L 2 (r_), and a = a(x) is a bounded measurable function 
on fl such that 0 < a 0 < a(x). While this is the simplest of hyperbolic conser- 
vation laws, solutions to (2.1) may contain discontinuities along characteristic 
lines x(s) defined by = f3. Solutions to (2.1) belong to the space of functions 
!/(fl) = {u G L 2 (fl) | vg G T 2 (fi)} where vp = (3 ■ Vu. 

2.2 Notation 

Throughout this work, notations and conventions standard in the literature 
on the mathematics and application of finite elements are used. Particularly, 
H S (Q) denotes the usual Sobolev space of functions with distributional deriva- 
tives of order s in Z, 2 (fi), equipped with the norm, 



Other notations and norms are defined in this section and where they first 
appear in the text. 

The starting point for the discontinuous Galerkin method is (2.1) defined 
on a partition of f l. Let Vh denote a partition of f) into N = N(Vh ) subdomains 
K with boundaries dK such that 

(i) N(V h ) < oo 

(ii) n = U {K : K G Vh] 

(iii) For any pair of elements K , L G Vh such that K ^ L, K fl L = 0 
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(iv) K are Lipschitzian domains with piecewise smooth boundaries. The 
outward unit norm to dK is denoted by n^. 


(v) dK _ = {x £ dK | (3 ■ uk < 0} and dK + = dK \ dK _ and no 
boundary dK coincides with a streamline, n/v • /3 / 0 

(vi) T* = U ^=1 dif D r_ coincides with T_ for every h > 0 


(vii) Tkl = dK D dL is an entire edge of both K and L 

(viii) The elements K 6 Vh are affine maps of a master element K = 
[—1, 1] x [—1, 1], K = F k (I<) as illustrated in Fig. 2.1. 


(ix) Vh € T where T is a family of quasi-uniform refinements. Let hx = 
diam(A" ) and px denote the supremum of all spheres contained in 
K; then for all Vh G T, there exist positive constants a and r, 
independent of h = max/f € p h hx, such that 

~t~ < r and — — < a (2.3) 

hx Pk 


We extend V(n) to the partition using the broken space 


I 


V{Vk) = n V{K) 

Ker h 

V(I<) = {veL 2 (K)\v p eL 2 {I<)} 


which admits discontinuities across element interfaces and use the following 
notations concerning functions v , w £ V(Vh) : 

v inlK = v| K (x), xedK 
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,ext K 


= y| L (x), X e dl< n dL 


v ± = limi;(x±e/3) 


(v,w) y = I vw\j3 • n 7 | ds, 7 C dK 


(i v ))l = («, v )i 


K w) K = f 

j Ti 


vwdx 


(2.4) 


MU = (»>«), 


We define the following norms for functions v £ V(Vh)' 


1 IMIU.„ = {iMli + IMli + «»»!*} 8 

{|MI 1 + IMI 1 + «» + ))L_ + «» _ )> 3 ff*} 


ll yfliK 


= {IMIJ + IMP* + <(« + - r_ + 


def 


E IIMII 

.Kef'll 


2 

BtK 


def 




def 


E 

,Kev h L 

£ 

,K&> h 


h,K 


Mil + ( 1 + i'ir I «» + )>L. 


Pk. 


MI1 + (i + ^)«0) 


2 

dK + 


) 


def 


E [«*IMI 1 + II»„* 


,Kev h 
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+ (1 + 6 ^-) («»+ - r. + 

where 0 is used to indicate the value of the parameter 

If 9 = 8hp then 0 K = 8 ^ 

Pk 

HO = hp then 6 r = 

‘ Pk 

If 6 = h then 9 K = h K 

If 6 = 0 then 0 K = 0 

If 0 = 1 then 0 K = 1 

The parameter p K appearing in the definition of d K in the mesh-dependent 
norm (2.5) will later represent the spectral order of the polynomial approxima- 
tion in K . The case in which the coefficient 9 h . = ^ in (2.5) plays an important 
role in the stability and error of the method, as we show later. Throughout C 
is used to denote a generic positive constant, not necessarily the same at each 
occurence. 



2.3 Weak Formulation 

Property (v) of the partition implies that solutions u £ V(Sl) to (2.1) are con- 
tinuous across element interfaces. Since the broken space admits discontinuities 
along element interfaces, we have the following problem corresponding to (2.1) 
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on the partition Vh- 

Find u 6 V(Vh) such that for every K £ Vh, 1 


up + au = f in K 

► 

u mt a (3 - UK = u ext h f3 • n k Vx £ dK _ \ dfl 

u int A /3 • n/c = </ /3 • rift- Vx £ dK- fl T_ 
for which weak solutions are sought satisfying 

Find u £ V’('Pfc) such that f° r every K £ Vh, \ 


(2.6) 


+ = (/,w) K WweV(K) 


(2.7) 


(u lnt ^ , v)dK_\9Q = {u ext h ,v)dK-\dQ Vu £ V(K) 

(u int K ,v) dK _nr- = { g,v) dK _nr _ Vu £ V{I<) 
where we have taken the absolute value of /3 • rift- for convenience. Next, we 
introduce a global parameter 6 which has a value of either 0 or 1. Recall that 
for any v £ V(K) we have that vp £ L 2 (K ) so that we can set w = v + S^vp in 
(2.7) and then add the boundary integral equations multiplied by any constant. 
It will be convenient to choose this constant to be (1 + and to write the 


method in the following abstract form. Let 


B k (u,v) = f (up + au,v + 6-^-vp) K + (l + 6-^-)(u + 

Pk Pk 


u 


, v + )dK-\r- 


+ (1 + f>^~)(u + y v + )dK- nr_ 

Pk 


( 2 . 8 ) 


L k {v) = {f ,v + 6^-vp) K + (l + 8^-){g, v) dK _nr- (2.9) 

Pk Pk 

where, by definition, u + = u int K and u~ = u ext K on dK-. Summing over all 
the elements in the partition yields the variational boundary value problem for 
weak solutions to (2.1) on the partition: 
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Find u £ V{Vh) such that 

B(u,v ) = L(v) , for every v £ V[Vh ) (2.10) 

where 


B(u,v) = £ B i<( u i v ) 

Kev h 

(2.11) 

L(v) d = f £ L k (v) 

Kev h 

(2.12) 


Remarks: 

(i) The case 8 = 0 in (2.8) is referred to as the ’’standard” discontinu- 
ous Galerkin method which can be viewed as a standard Galerkin 
method for a single element with weakly imposed boundary condi- 
tions for elements lying on the inflow boundary and weakly imposed 
continuity for elements on the interior of the domain. 

(ii) The case 8 = 1 in (2.8) is the hp extension of the so-called ’’stream- 

line upwind” discontinuous Galerkin method [29]. The modification 
of the test function is important when approximating solutions with 
sharp gradients as the additional term in the test function adds 
diffusion in the streamline direction without modifying the conser- 
vation law, i.e., without destroying accuracy in regions where the 
solution is smooth. I 

Lemma 2 Let the bilinear form B ( •,•) be defined by (2.11) and (2.8). Then 
there exists positive constants a, M\, and M 2 , independent ofh K andp K , such 
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that 


B(v,v) > a|||u|||^ p/ j 

B(w,v ) < Mi ( J2 (l + 6%)| 

l Kzv h \ Pk ) 


(2.13) 


w 


n,p,K 


X < V 


h k 


+ IMIi + «* + )> 

Kev h PK 


+ \\2 

dh\ 


(2.14) 


«(»,») < Ml E IIWII1 + *E ^[IKIII + ((»■)> 


Kev h 


K€P h PK 


— \ \ 2 
3A + 


X 


£[1+4 


.A'6P„ 




%.») < MA £ ( h -<^) iiwii 2 „ 

[h ev^K Pk J 


) 


(Kev h 

for every w,v £ V(Vh)- 


Pk 


(2.15) 


(2.16) 


Proof: (i) From the definition of B ( •, •), 

S(n,n) > min(l,niina(x)) £ ( (l + <5^-1 (p,p/j) a -+ + IM 


+ (l + <5% | [{{v 
V Pk 


Kev h 

+\\ 2 


Pk 


Pk 


))aK_\r_ ~ ( v , v+ )dK-\ r_ + <( u ))L'_nr_) | 

Equation (2.13) follows by substituting the results of applying Green’s formula 
to the term ( v,vp ) K , that is, 

E (». ”/>)* = ;E («®'»5jf-\r- + «o + ))L_\r_ + «v)>lirnwi) (2-17) 


*€P h 


‘ft'ePh 
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into the above and choosing a = | min(l, min xe n a(x)). 

(ii) Applying the Schwarz inequality to £?(•,•) as defined in (2.11) and (2.8) 
yields 

B(w,v) < iMkn £ IWL-lklU + ^IMLIklL + IHUIIHI, 

Kev h L Pk 

+ ^IMUML + (l + ^pj^) + ~ W ))dK_\r_({ y+ )) 8 /»_\r_ 

+ ( 1 +*^r)«“’ + )) 8 K .nr_«» + » 3 K -nr. 

< II-IU.J E fi + 40 [iMli + IK- 

[ Ker h \ Pk / 

+ ((« ,+ — w ))a/f_\r_ + ((^))aft'nan] } 

x III 2S ^t\\ v p\\1 + 2 IMIk + f 1 + (( v+ ))Ik. } 

U'eP h L Pk \ Pk J JJ 

Equation (2.14) follows by selecting M\ — \/2||«||oo,n- 

(iii) Equation (2. 15) is obtained by applying Greens formula to the term (re, vp) K 
and (w/ 3 , v) K in (2.11), and applying the Cauchy-Schwarz inequality. 

(iv) Equation (2.16) is obtained by adding (2.14) and (2.15), applying the 

Cauchy-Schwarz inequality to the result, and selecting M 2 — |||a||oo,fi- 1 

Corollary 1 Let 8 = 1 in (2.8). Then 

B(v,v) > a\\\v\\\l p0 (2.18) 

and there exists a constant r 0 such that if < r 0 V/f € Vh, 


B{w,v) < A/ 1 '|||u ;||| 1i/ 3 


lkv h PK 


1 

2 


(2.19) 
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B(w,v) < M , 1 \\\\w\Wl+ E ¥|KIIk1 IIMII (2-20) 

l Kev h PK ) 

B(w,v) < M 2 (l+r 0 )\\\w\\\ B \\\v\\\ B (2.21) 

Proof: Set 6 = 1 in (2.13)-(2.16) and choose M[ = Mi\/1 + romax(l, yjr o). I 


Corollary 2 Let 6 = 0 in (2.8). Then 


B( v, 

v) 

> 

ollMII^ 

(2.22) 

B{w , 

v) 

< 


(2.23) 

B{w, 

v ) 

< 

M.IIMIUIIIHIII.jj 

(2.24) 

B{w , 

v ) 

< 

m 2 ||WIUIHII» 

(2.25) 


I 


Remark: 

Note that modified test function for the streamline upwind discontinuous Galerkin 
method results in improved stability of the bilinear form when compared to the 
standard Galerkin method (see (2.18) and (2.22)). The coercivity of the bilin- 
ear form for the streamline upwind discontinuous Galerkin method contains 
||v/9||#c terms which do not appear in the coercivity condition for the standard 
Galerkin method. The significance of this additional stability is less important 
as approaches zero. ® 


2.4 hp Finite Element Approximation 

We seek approximate solutions to (2.10) in the finite dimensional subspace 
V p (Vh ) C V(Vh) defined as follows: 


26 


V,(Pk) = {» £ 1) : v\ K = «|„ o F„ € Q’ K (K)) (2.26) 

where Q A ( K ) is the space of tensor products of polynomials of degree p K 
defined on the master element. The basis for Q Pk (K) is formed by tensor 
products of Legendre polynomials. We use the notation v K £ Q P,< (I{) to mean 
v K £ Q Ph (K) and v K = v K o F K . We have the following inverse estimates for 
polynomials on a single element: 

Lemma 3 Let I( £ 1Z 2 be an affine map of a master element K = [—1,1] x 
[—1,1]; that is K = Fk(I\). Let 7 denote any edge of dK which is an affine 
map of a master edge 7 = [—1,1], Let Wr be a polynomial of degree p^ defined 
on the master element. Let wk — o Fk denote the image of wk under the 
transformation Fk- Then (3 ■ 'Vwk satisfies the following: 

W& VwkWk < cf£|K.|| A . ( 2 . 27 ) 

tlx 

W-Vv, K )) 1 < C V - (2.28) 

n K 

where the constants C are independent ofh K ,p K , andwx- 

Proof: For polynomials of degree p K on the master element (see Dorr [21]): 

5: Cp 2 K 1 1 U>K 1 1 K 
I^A'U-y 5: < CpJrll^A'H'Y 


where the constants C > 0 depends on s, but not on p K or Wr. 


(2.29) 

(2.30) 
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For affine mappings Fk, a standard scaling argument (see Ciarlet [11]) 
yields that for an integer s > 0, there exist constants C > 0 such that 


Wk\s,K < 


(2.31) 


Ch^ s |u>a| s , 7 

(2.32) 


ChZ'lwK :U 

(2.33) 

|^A'|s,7 — 

Ch s K 2 |iwa-| S i7 

(2.34) 


where C depends on s, a, and r (see (2.3)), but not on h K , p K , or w k- 

The first estimate (2.27) follows by combining (2.31), (2.29), and (2.33). 
The second estimate (2.28) follows from (2.32), (2.30), and (2.34). I 

Lemma 4 (Babuska and Suri [3]) Let I< G Vh, 7 denote any edge of dK, 
and u € H S (K). Then there exists a constant C — C(s,r,a) independent of 
u,pk, and h K , and a sequence z v h € Q” K (K),p K = 1,2,..., such that for every 
0 < r < p K , 

ll«-^IU < S>0 (2.35) 

Pk 

ll«-*sik, < C^\\u\\ StK , s>l (2.36) 

Pk 2 

where u = min(p^' + 1, $). I 

The approximate solution to (2.10) is obtained by replacing the exact 
solution u G V(Vh ) by u p h G V p (Vh ) and the test function v G V{fPh) by 
< € V p (V k ): 


Find u p h G V p (Vh) such that 
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fl«, »I) = i«) , V»J € V p (V h ) (2.37) 

The improved stability of the streamline upwind discontinuous Galerkin 
method, 6 = 1 in (2.8), is recovered by the standard discontinuous Galerkin 
method, 8 = 0 in ( 2 . 8 ), on the finite dimensional space V p (Vh,). 


Lemma 5 Let 8 = 0 in (2.8). Then for every v v h € V p {Vh) there exists a 
w h ^ V p (Vh) such that 

B(vlwl)>a'\\\vl\\\l p<p (2.38) 

and 

IIKIIUwi < C|IKIII»p.O (2-39) 

where the positive constants a' and C are independent of h K , p K , and v p h . 


Proof: Define the restriction of w v h G V p (Vh) to an element K G Vh as 

<\ K = <\„ + 7^/3 • V<L- (2.40) 

Pk 

where 7 6 (0, 1] is defined later in the proof. Dropping the /i,p, and I\ scripts 
for ease in notation, we have 

B k {v,w) = f (vp + av)(v + J- 5 -vp) dx 
Jk Px 

+ f (v + - v~)(v + + ')^-v (} + )\(3 ■ n K \ ds 
JdK.\ r_ p K 

+ f v + (v + + i^-vp + )\f3 • n K \ ds 

JdK-n r_ px 

a o|MI* + 7%IMI* + [ vvp dx 
Pk Jk 


> 
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+ 1 pj^ J K VV P dx + {( t ’ + ))9A-_nr_ + (( v+ ))dK-\ r_ 

- [ v + v~\/3 ■ n A -| ds + 7 % [ (v + -v~)v j3 ' i 

Jdh'A r_ p 2 K JdK-\ r_ 

+ 7 ~T I v + v p + \f3 • rift- 1 ds 
Pft 7aA'_nr_ 

where a 0 = min xe n «(x). Noting that 

L VV P dx = \ I (0 2 |/3 • haI ds-)- ( (u + ) 2 |/3 • rift- 1 ds 

vft ^ JdK + ^ JdK- 

and that from Lemma 3 

| j W /3 dx\ < Cl ~^\\v\\ 2 K 

JK tlx 

I J dK _ v+v P + \fi • nft'| ds | < c 2 ^-({v + )) 2 dK _ 

we have 

B k (v,w) > (g 0 — c i T ) 1 1 ^ 1 1 A' + l~^~\\ v p\\]< 

Pk 

+ 2 (( u+ ))aA'_\r_ + ~ C 2 7)(( u ))aft-_nr_ 

- f v + u _ |/3 • nft- 1 ds + 7 ^ / (w + - v - )^ 4 

JdK-\r~ p 2 K JdK-\ r_ 

Using the Schwarz inequality and the previous inequalities, one can 

/ K -\r “ i; “)^ + l/ 3 'nA'| *| < ^7 (((v + ))L_ 

+ (( v ))aA'_\r_) 

Now summing over all the elements K € Vh and realizing that 
£ <5«0>L_\r_ - r. 


1/3 ■ rift- 1 ds 
(2.41) 


|/3 • nft- 1 ds 
show that 

\r_ 
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+ (^ “ ^7)<(v + ))9 A -_\r_ - J dKAr _ v+v I/ 3 ' n *'l <M 

> |((®)>? + +min(l i-^y) £ <( u+ “ 0>L'_\r_ 

z z z Kev h 

results in 

«(»;,<) > (*k.-c, 7 )|HIJ + 7 £ tHMIk 

Kev h Pk 

+ (j-C27)«»)>?_ + j(W>? t 

+ min(l ,5-^7) £ (( v * - v ~))dK-\r- 
z Kev h 

Choosing 7 = min(i, g^-) yields the first inequality. 

The second inequality easily follows from the definition of and Lemma 3. 

I 


2.5 A Priori Error Estimate 

The discontinuous Galerkin method (2.37) was first analyzed by Lesaint and 
Raviart [33] for a given fixed value of p K , i.e. for the case in which p K = p for 
every element K 6 The error in a solution Uh to (2.37) approximating an 
exact solution u € H s (ft) to (2.10) was shown to be 

\\u-u h \\n<Ch*- x \\u\\ afl 

This estimate is not optimal in the sense of interpolation error estimates and 
was improved by Johnson and Pitkaranta [30]. Using a mesh-dependent norm, 
they showed that 


\\\u-u h \\\ h , 0 <Ch°-t\\u\\ s , Q 
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While this estimate is not optimal in the sense of interpolation error estimates 
for ||e||n = ||u — u/i||n, it is optimal with respect to J~h^\\ep\\K and {{e + — 
e_ ))aA'_\r_- We shall derive estimates similar to Johnson and Pitkaranta [30] 
taking into account that p K is not constant. 


Theorem 1 Let u € H s (fl) be a solution to (2.10) and let u p h be a solution to 
(2.37). Then there exists a positive constant C, independent of h K , p K> and u, 
such that the error, e = u — u p h) satisfies the following estimate 




(2.42) 


where u K = min(p A . + l,s). 


Proof: Let ll v h u G V p (Vh ) be an approximation of u that satisfies the estimates 
in Lemma 4 and write 

e = u — u p h = u — n£u + Il£u — u p h (2.43) 

which implies that 

lll e IIUp,/? — lll«-nM||;^ + ||K-nMiu P ,, 

= f IIMIUp./? + IIMIUp,/? (2.44) 

where, to simplify the notation, we set r/ = u — Il£u and w = u p — H p h u. 
Realizing that 

{Kev h Pk \ PkJ J J 
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combined with Lemma 4 yields bounds for the first term in (2.44): 


\h P ,0 < C) E 

(Kev h 


2 tl Tf 1 

Ik 


1 


j2^ k .2ii k 

, 11 K , r f*K 

n 2 s ' d 2 s -1 ' 

” K " K ” K 


Nil 


< Cl E >T-ax 
\Kev h P K 


1 hj£ 

? 2 
Pk Pff 


(2.45) 


where fi K = min(p K + 1,5). Bounds for the second term in (2.44) follow from 
the orthogonality condition which is obtained by subtracting (2.37) from (2.10): 


B{e,v p h ) = B{r],v p h ) - B(w,v p h ) = 0 , Vv%.€ V p (P h ) 


We choose v v h = v 8 in (2.46) where the particular choice for v 5 depends 
on the parameter 5 in (2.8). For S = 0, we choose v s to be the function which 
satisfies Lemma 5. For 6=1, we choose v 6 = w and combine (2.46) with 
Corollary 1. The result for either case is 

C\\W\\l r ,e<B(wy) = B( V y) (2.47) 

Integrating the terms (f]p, vS ) K and (Vi v p) K by parts in the definition of f?(-,-) 
in (2.11) yields 

< IMkn E [mIkII”‘II« + ^MIk-^KIIk 

Ker h l \Jh K Pk 

Pk Pk 

+ ))9A'+\r + (( y5+ — v ))aA' + \r + 

+ ))dK+nr + (( vS ))dA' + nr + 

< ||«IUn(l + <5)( E \( l + r)^w 2 K 

{Kev h i\ a K) 
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(2.48) 


Recall that for our choice of V s we have |||v 5 ||Up,/ 3 = |||Hlk P ,/? when 5 = 1 
and 1 1 1^| | \hp,/3 < CHHIliW when 5 = 0. Equations (2.48), (2.47), and the 
estimates in Lemma 4 imply that 


HI \h P ,p < 


h]- K ~ l n 


+ 


< 



+ 5 1+5 


hK 

Pk. 


. c h 2 ' 

2 , 1 ,—, 

Pk P k PI ) 


« 


5 1 A' 


(2.49) 


Combining (2.49), (2.45), and (2.44) completes the proof. 


I 


Remarks: 


(i) For < 1, the estimate becomes |||e|||^ Pi; g < C{£,K&> h ^Sf =a r IMI?,y}* 

(ii) For p K =constant, the a priori error estimate reduces to the one 
derived by Johnson and Pitkaranta[30]. 

(iii) Let h = max KeVh h K and p = min p K , then |||e||| fcPi/ , < ^r||u||,,n. 

I 


2.6 Implementation Issues 

In the preceeding sections of this Chapter, the discontinuous Galerkin methods 
were represented as global methods for the purpose of analysis. The approxi- 
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mate problem is actually a local one since the approximate solution in an ele- 
ment is independent of the solution in the neighboring elements with the only 
coupling between elements occuring weakly through the fluxes on the element 
inflow boundary. Assume that u p ~ is known on dK _, then the approximate 
solution in element K satisfies 




) = i,KL),Vi€«"(A') 


h\K^ u h\K 


(2.50) 


where 


KKA) = <„),< 


Pk 


+ (i + *^)K + ,< + W- 

Pk 


(2.51) 

(2.52) 


L K ( v h) — {f, v h + 8-^-v p h0 ) K + (1 + 8-^-)(u v h ,v p b + ) 9K _\ r _ 

Pk Pk 

+ (1 + 8-^-){g,v p h J *') dK _ ri Y- (2.53) 

Pk 

In order to solve (2.37) in this fashion, one must define an ordering of elements 
that starts at the domain inflow boundary and sweeps through the partition in 
such a way that u p h is known on <9A_ prior to solving (2.50). Such an ordering 
always exists (see [33]) and is fairly straightforward to construct. This is the 
optimal solution technique for solving the linear model problem where element 
inflow boundaries can be identified a priori. However, an alternate approach is 
needed for solving the conservation law (1.5) where the fluxes depend on the 
solution, and thus, element inflow boundaries cannot be identified a priori. 


With the aim of solving more general problems in mind, we will solve the 
linear model problem in a way that is easily extendable to the nonlinear case, 
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that takes full advantage of the discontinuous approximation, and is amenable 
to parallel computations; that is, by solving the time-dependent conservation 
law for the steady-state solution. Since time accuracy is not important in 
obtaining the steady solution, we use the classical forward or backward Euler 
time marching with a truncation error of 0(At 2 ). Let u n+1 = u p h (-,t n+1 ) where 
t n + 1 = (n + 1)A t and A t is the time step increment. Assuming the solution at 
time t n is known, then the forward Euler version of the scheme is given by 

£ (»” + ',») k . = £ («“,»)« + Ai[L(»)- B(«»l (2.5-1) 

!<e'P h K<zV h 

and the backward Euler version is given by 

£ K + \v) a . + A t £ BJu"*\v)= £ („",„)* + A t £ l K (v) (2.55) 

Kt-Ph KeVh Kev h Kev h 

To preserve the local character of the method, the inflow boundary terms ap- 
pearing in the definition of (see 2.53) are evaluated at time level t n . 

The initial data, u° = u^(-,0), needed to complete the initial-boundary-value 
problem is taken to be a uniform field with a value associated with the inflow 
boundary conditions. 

A time-accurate Runge-Kutta time marching scheme for nonlinear con- 
servation laws is described in detail in Chapter 7. It can essentially be written 
as a sequence of steps in the same form as (2.54). Parallel implementation of 
the time marching discontinuous Galerkin methods for general hp meshes is 
described in Chapter 6. 



Chapter 3 


A Posteriori Error Estimation 


The a priori estimates derived in the previous chapter are useful for predicting 
how the error in numerical solutions behaves with h-refinement or p-enrichment. 
Unfortunately, their usefulness in assessing the accuracy of a given numerical 
solution is limited since the estimate involves unknown constants and the exact 
solution we are approximating. Nevertheless, a priori error estimates such as 
(2.42) and interpolation error estimates such as (2.45) have been used exten- 
sively as error indicators to drive adaptive methods for hyperbolic problems 
[16], [43], [36], Typically the unknown constant is set to unity and some post- 
processing of the approximate solution is used in place of the exact solution. 
While the element contributions to these global estimates may provide some 
relative measure of the local error, this approach in general fails to provide a 
reliable estimate of the actual error in a particular numerical solution and can 
be grossly in error. 

In this chapter, we derive error estimates which are computed locally 
on a single element and contribute to a global error estimate which is accu- 
rate enough to provide a reliable assessment of the quality of the approximate 
solution. 
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3.1 Element Residual Method 

The estimates derived here, based on the element residual method, are simil- 
iar to those proposed by Bank and Weiser [4] for elliptic problems and Oden, 
Demkowicz, Strouboulis, and Devloo [38] for solid and fluid mechanics prob- 
lems. The element residual method was extended to /ip-approximations for 
elliptic problems by Oden, Demkowicz, Rachowicz, and Westerman [37]. A 
global estimate of the error is obtained by summing element indicators which 
are the solutions to a local problem with the element residual as data. In 
references [38] and [37] , the local problem is of the same form as the global 
problem. 

For continuous finite element approximations, the element residual in- 
volves fluxes on the boundary of an element. Since the fluxes are multi-valued, 
an averaged flux is used. Recently, Ainsworth and Oden [1], [2] have shown 
that it is possible to use a self-equilibrating average flux that results in an error 
estimate which is equivalent to the actual error and can be asymptotically ex- 
act for certain elliptic problems. For a discontinuous approximation, the jump 
in the element boundary flux arises naturally in the residual, eliminating the 
need for flux balancing. 

The main difficulty with our formulation for hyperbolic conservation 
laws is that the norms associated with the continuity and coercivity of the 
bilinear form are different. Therefore, use of different norms makes it impossible 
to construct a single local problem which results in an upper and lower bound 
of the error in the same (or an equivalent) norm. 

In the following sections we show that it is possible, however, to con- 
struct one local problem with a solution that provides a lower bound on the 
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actual error and another local problem with a solution that provides an upper 
bound on the actual error. We also show that a local problem based on the 
original problem results in a local lower bound. Moreover, if the approximation 
of the solution to this local problem is limited to a certain class, then the es- 
timate is equivalent to another commonly used approach: estimating the error 
as the difference between a newly contracted (and hopefully more accurate) 
solution and the approximate solution on hand. 

3.2 A Global Lower Bound on the Error 

We define a local problem which results in a lower bound on the error in a 
sense to be defined precisely later. Let u K € Q A (I\) denote the approximate 
solution in element K and ip K 6 V(K) be the solution to the following local 
problem, 

= B K (e K ,v K ) = L k (v k ) - B k (u k ,v k ) Vv k e V(I<) (3.1) 

where 

K (Pk> v k) = ^r(0 ■ *¥>k>0 ■ Vv k)k + a(<p K ,v K ) K + (<p k ,v k ) 9K _ (3.2) 

Pk 

and a > 0 is a constant. Then 

d I (^,v)= £ A L k {v k ,v k ) (3.3) 

K€V h 

induces a norm on V(Vh) which will be referred to as the /l £ -norm: 

M\ =A L (v,v) Vu E V{V h ) (3.4) 

The solution to the local problem (3.1) provides a lower bound on the 
error in the following sense: 
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Lemma 6 Let p> £ V{Vh) be the solution to the following problem: 


A {<p,v) = B(e,v ) Vu £ V(Vh) 


(3.5) 


There exist positive constants k\ and r 0 st/c/i that if < ro V/L € Vh then 

P/v 


IMI.L <*iIIHIIij» 


(3.6) 


Proof: |bl| 2 , = 4 (v?,<p) = J5(c,v>) 


< Mi £ t+« 




Pk 




X 


V 

K<TPh L Pk- 


+ IK- 




dK- 




om (2.14) 


< M x (l + fr 0 )max(l,- 7 =)|||e||| ||v?| 


The desired inequality (3.6) follows by choosing Aq = Mi(l + <5r 0 )max(l, As). 


3.3 A Global Upper Bound on the Error 

For simplicity, the estimates in this section are derived for the case when 8 = 1 
in (2.8). We construct a local problem which results in an upper bound on the 
error. Let u K £ Q* K (K) denote the approximate solution in element K and 
tp K be the solution to the following local problem, 


A kWk> v k) = B A e K> v tc) = L k( v k) - B k{u k ,v k ) VveV{I<) (3.7) 
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where 

= %r(0 ■ ■ Vv,), + a(^,v K ) K (3.8) 

Pk 

and a > 0 is a constant. Then the /l^-norm is defined as 

IMI% = = IC KMk^k) (3-9) 

K£P h 

The solution to the local problem (3.7) provides an upper bound on the error 
in the following sense: 

Lemma 7 Let ip 6 V(Vk) be the solution to the following problem : 

A U (ip,v) = B(e,v) VveV(Vh) (3.10) 

where 8 = 1 in the definition of 5(*,-) in (2.11). Then there exists a positive 
constant &2 such that 

\M>\\ a u >k2\\\e\\\ hPt( ] (3.11) 

Proof: Using (2.18) of Corollary 1, 

a lll e llliU® - B ( e ’ e ) = ^ £/ (V’> e ) 

< IWUMI „ 

A A 

kev^Pk 

< max(l, V5)||^|| { Y. IMIJ + INI*)* 

Kev h Pk 

< max(l,v^)||^|| „|||e||| fcPt/J 

A 

Choosing &2 = — f /rr. completes the proof. I 

rn&xi x j y cl j 
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3.4 A Local Lower Bound on the Error 

Recall the bilinear form (2.50) which characterizes the space marching form of 
the discontinuous Galerkin method: 

Pk 

+ (i +s^)K\vf) sh -_ 

Pk 

Introducing a local norm, 

IKH., = felKX + IMli 

h L Pk 

/ \ i 

+ j( 1 + «^|)((''“»5 Kt ’ (3.12) 

and using Green’s formula, it is easy to show that there exist a constant C > 0 
such that 

KK^ k ) > C\\v\\l Vv K 6 V(K) (3.13) 

a K 

Now consider the following local problem: 

Find ip K £ V{K) such that 

= Bk^k^k) . € V{I<) (3.14) 

then <p K provides a local lower bound on the error in the following sense. 

Lemma 8 Let y> K E V (K) be the solution to ( 3.1 4)- Then there exists a 
constant k 3 > 0 such that 

\\¥k\\b k - K + s ^j\\\ e \\\i,P,K (3.15) 
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Proof: Setting v K = tp K in (3.14) and using (3.13) yields 


C\Wk\\\ < B K (<p K ,<f K ) = B K (e K ,<p K ) (3.16) 

Setting M k = max(l, IMIoo,*) an d applying Young’s inequality, ab < £ a 2 + 
eft 2 , e > 0, to each term in B h .(e K , tp K ) yields 


B k{ c k^k) < 


Mk 

4e 


.h 


l+^j\M\l e , K + 2M K e\\ VK \\] 


(3.17) 


Selecting e < in (3.17) and combining with (3.16) completes the proof. I 


3.5 Approximation of the Local Problems 

An approximate solution to the local problem measured in the corresponding 
norm serves as a local error indicator for an element. Since the discontinuous 
Galerkin solution satisfies the orthogonality condition, 

B K (e,v) = 0 Vu € Q PK {I<) (3.18) 

we must approximate the error indicator with a polynomial of degree p K + a K 
where a K > l in order for the discrete local problem to have a non-trivial 
solution. If a complete polynomial of degree p K + a K (on the master element) 
is used to approximate the solution to the local problem, then the discrete local 
problem requires the solution of a system of order ( p K +cr K + l) 2 . This system 
can be fairly large compared to the system of order ( p K + l) 2 equations used 
to obtain the approximate solution for which we are estimating the error. 

Since ( p K + 1) 2 terms on the right hand side of the discrete local problem 
(corresponding to (3.18)) are zero, we can make a simplification by approxi- 
mating the solution to the local problem in the space Q p x +tT K (K ) \ Q p x(K). 
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In other words, the solution to the local problem can be approximated with 
incomplete polynomials of degree p K + a K by neglecting the terms associ- 
ated with polynomials of degree p K . This simplification results in a system 
of cr K (a K + 2 p K + 2) equations for each element. 

The size of the local problem can be further reduced by approximating 
its solution using only the ’’bubble” functions in the enriched space denoted 
by Qq k ^ <Tk (K) \ Qq k (I\). These are the polynomials in Q p i<+ a K (/i) \ Q p i<(K ) 
which are zero on the boundary of an element. This additional simplification 
results in a system of cr K (a I<; + 2 p K — 2) equations which is smaller than system 
of equations used to obtain the approximate solution. 


3.6 Remarks Concerning an Alternate Approach 

Suppose that an approximation U to the exact solution u can be constructed 
which is more accurate than the approximate solution on hand, u p h . Then a 
simple estimate of the error, e — n — u v h ~ u ~ U + U — is 9 = \\U — u£|| 
where || • || is any suitable norm. Using the triangle inequality, we have 


||e||-||u-l7|| <0<||e|| + ||u-t/| 


or equivalently 


, < JL< I, IN -i'll 

I I I I — il ll — ~T* II ii 


ll e ll ll e ll ll e ll 

If ||u — U || << | |e 1 1 = ||u — <||, then 6 = \\U — <|| is a good error estimate 
with an effectivity index near unity. The main difficulty with this approach 
is to efficiently construct such a U. One obvious strategy for constructing a 
more accurate approximation, U, is to re-solve the approximate boundary-value 
problem on an enriched space. 
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For continuous finite element approximations, this leads to a global sys- 
tem of equations that is much larger, and therefore much more costly to solve, 
than the original problem. For a discontinuous approximation, re-solving the 
problem on an enriched space of complete polynomials is still more costly than 
the original problem, but is no more costly than solving the local problems in 
section 3. 2-3. 4 on the complete polynomial space. 

The computational cost of this approach can be further reduced by 
’’freezing” the lower-order solution and re-solving the problem on an incomplete 
polynomial space of bubble functions. In other words, let U K — u p h \ K +w K where 
w K € Q P 0 k+<Tk {K) \ Ql K {K) satisfies 

B k (w k ,v) = L k (v) - B h .(u p h ,v) \/v e Ql K+a,< {K) \ Q P 0 «(K) (3.19) 

In this case, the error estimate is 0 = \\U — w£|| = ||t»||. Note that this is equiv- 
alent to solving the local problem (3.14) on the space of bubble functions. We 
remark that Peraire and Morgan [44] simply post-process the approximate solu- 
tion to obtain the degrees-of-freedom (higher-order derivatives) corresponding 
to the bubble functions. 



Chapter 4 


An hp - Adaptive Strategy 


The /ip-adaptive strategy used here is based on a 3-step strategy developed by 
Oden, Patra, and Feng [39]. The strategy was developed for a large class of 
elliptic problems and has been shown to yield exponential rates of convergence 
with respect to CPU time [39]. The hp - adaptive strategy is based on a reliable a 
posteriori error estimate for determining the error in the approximate solution 
and an a priori error estimate for determining how to modify the mesh to 
improve the solution accuracy to a specified level. The goal of the /ip-strategy 
is to deliver a solution with a specified error in only three steps: 

(i) Construct an initial partition Vo containing N(Vo) elements. The 
elements in Vq can be of uniform p K = p 0 and essentially uniform 
in h K = h 0 . Solve the problem of interest on V Pq (Vq) and estimate 
the error. 

(ii) Construct a partition V\ by subdividing each element in Vo into 
the number of elements required to equally distribute the error and 
reduce it to a specified level. Solve the problem on V p (V\) and 
estimate the error. 
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(iii) Enrich the approximation space by increasing p K for every K £ V\ 
in such a way to equally distribute the error in smooth regions and 
reduce it to the specified level. Solve the problem on the enriched 
space V p i (Vi) and estimate the error. 

If the estimated error in the solution is larger than the specified error after the 
third step, then it is necessary to repeat steps (ii) and (iii) until the desired error 
is attained. For discontinuous solutions, p-enrichments in step (iii) are confined 
to elements in regions where the solution is smooth, since higher-order elements 
at discontinuities may result in oscillatory solutions. Moreover, p-enrichment 
of elements in regions where the solution is of low regularity does not improve 
the accuracy of the approximation, as indicated by the riori estimate (2.42). 

The data structure used for the resulting ftp-meshes is based on the 
work of Demkowicz, Oden, Rachowicz, and Hardy[19] for continuous finite ele- 
ment approximations. The data structure for the initial mesh consists of nodal 
coordinates, element connectivities, boundary conditions, a list of neighboring 
elements, and element orders based on the number of edge and interior degrees 
of freedom. Refinements are achieved via bisection of an element in the initial 
mesh and are added using a tree data structure. The data structure routines 
enforce a mesh irregularity index of 1 and enrich the order of an edge for an 
element with a neighboring element of higher-order. These two properties are 
necessary for maintaining continuity of the finite element approximation and 
are not needed when using the discontinuous Galerkin approximation. The 
data structure for the degrees of freedom for the discontinuous approximation 
consists of three arrays to store the vertex, edge, and interior degrees of free- 
dom. For non-uniform p-meshes, the edge (and interior) degrees of freedom are 
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stored consecutively with integer arrays which provide the address of the first 
location of the edge (and interior) degrees of freedom in the global array. 

We make some basic assumptions before describing the adaptive strat- 
egy in detail. First, recall that the norm used in the a priori estimate (2.42) 
includes jump terms on the element inflow boundaries: 


E 


Kev h IPk 


+ ll< + (( e+ ~ e ))Ik-\ r_ + {( e ))lh-ndn 


< C 


E 

Kev h 


A h K 

max [ 1, — 


rfivK 
r k 


Pk 


WWh 


(4.1) 


For the a posteriori estimate we have several choices: the solution to 
the local problem which yields a lower bound (3.1), the solution to the local 
problem which yields an upper bound (3.7), and the simple estimate obtained 
by re-solving the problem on an enriched polynomial space (3.14). Since the 
norm induced by the lower bound local problem contains the error on the 
element inflow boundary, and not the jump in the error as in (2.42), we cannot 
use this estimate to drive our adaptive strategy. Fortunately, the norm induced 
by the upper bound local problem contains no element boundary terms and the 
simple estimate can be measured in any norm desired. We will use the norm 
induced by upper bound local problem with a = 1 and assume that the meshes 
at each step in the adaptive procedure are such that max ^1, = 1. Noting 

that (2.42) is valid if we drop the jump terms we can write 


' C ''au = 



£tIM 

Pk 


1 + IWI 


< c 



,2fi K 
l K 


2u 

Pis 


i-lK 


K 


(4.2) 
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We assume that the a posteriori error estimate is a good enough approxima- 
tion to the actual error to replace the left hand side of (2.42) and treat the 
inequality as an equality. Admittedly, this may not be a good assumption for 
coarse meshes and rough solutions, but it provides some means of predicting 
the structure of the new mesh. 


Since the adaptive procedure is based on refinement and then enrich- 
ment of an initial mesh, the term remains constant throughout the 

adaptive process and can be calculated from the estimated error on the initial 
mesh. Let 0o denote the estimated error for the solution uq 6 ^(To) s ^ e P 
(i). Then, 


el 


X el'K 

Ker 0 


2 »I< 
K 


= X ^c\\ 


U\ 


Ker o Po 


a 2 * 1 * 

2 v- n K 

r,K 2 v K lx K 

KeTo Po 


(4.3) 


or at the element level 


K 


El 

h>' 


700, 


K 


(4.4) 


4.1 The h - Refinement Step 

Let r] T denote the target error to be achieved by the entire 3-step strategy. We 
specify r] T as a percentage of the solution measured in the same norm as the 
error to assign it some physical relevance. For the A-step, we seek a partition 
which will deliver an intermediate target error r) h = ccr) T where a is some 
constant chosen so that t] T < rj h . The partition V\ is constructed by subdividing 
each element K € Vo into n K elements such that the error = 77 ^ ( 1 1« 0 1 1 + 0 O ) 

is distributed equally among the elements in Following (4.3), we have 

Tl WJ- Ti TS 1 2/i -wr 

*?= E E«h = E E%7 a 1 

KePo L = 1 Kev o L = 1 Po 


(4.5) 
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where 



K = E A l (4.6) 

L = I 


For a mesh which achieves an equally distributed error, 




/to 


, L = N(P,) 


(4.7) 


where 


N(V 0 ) 

N i'Pl) = J2 n K 

K = 1 

Combining (4. 5), (4. 6), and (4.7) yields 

= ^r^(Pi); /f = i,...jv(7>„) 


(4.8) 


(4.9) 


Equations (4.8) and (4.9) are solved iteratively to determine n K for each element 

I< € Vo. 

Remarks: 


(i) A value of n K < 0 signals that de- refinement is needed to equi- 
distribute the error. Although not implemented in this work, de- 
refinement significantly decreases the computational cost of the 
overall process. 

(ii) The largest local errors will occur in elements which contain a dis- 
continuity. These elements will receive the highest level of refine- 


ment. 
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(iii) For ease in implementation, the refinement of an element K G Vo 
is limited to 2 levels. 

(iv) The parameters p K = min(po + — |) and v K = r — 1 are global 

constants dictated by the regularity of the exact solution u G V(Q). 
Formally, (2.42) is not valid when the solution contains a discontinu- 
ity on the interior of an element. However, numerical experiments 
suggest that (2.42) the rate of convergence of the error is p K = | 
when the solution contains a discontinuity which is not alligned 
with the element interfaces. This value of p K is consistent with the 
finite difference results of Sanders [47]. fl 


4.2 The p-Enrichment Step 


Let 0 h denote the estimated error in the solution u 1 G V p (V\) obtained in step 
(ii). Treating the a priori estimate as an equality, 

K= E t* = E (4.io) 

Kev i KeVi Po 

which gives the constants A K onV\. 


t\„ = 


p 0 

*;*• 




(4.11) 


The error is reduced by constructing a distribution of polynomial orders, 
Pi, where the polynomial order of each element in Vi is selected to equally 
distribute the target error. Setting 0? = f?:r(||fii|| + Oh) and using the a 

A U 7 

priori estimate, we have 


0 2 t 


= E 

KeVi 
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V T,K 


i 2 ^ 
l K 

2u k 
K eVx Pk 


= E 


(4.12) 


51 


where for an equally distributed error 


Ot,k = 


yjN(T>y) 


(4.13) 


Combining (4.11), (4.12), and (4.13) yields the new value p K for each element 
in V\. 

„ Po • t h,Ks/N(K) 

ry K — 1 


PS = 


(4.14) 


For smooth solutions, the parameter p K in (4.12) depends on p K which 
is unknown at this point. Here we use the value which is actually associated 
with po- More discussion on the parameters p K and nu K is given in section 4.4. 


For discontinuous solutions, p-enrichments are confined to regions where 
the solution is smooth since increasing p at discontinuities may result in os- 
cillatory solutions and does not improve the accuracy of the approximation. 
The local regularity of the solution on each element K € Vo is estimated by 
computing the rate of convergence, fi K , of the local error in steps (i) and (ii): 

. = lQ g0°,ft' ~ Mg ^2=1 °l,L 
^ log h K - log 

From the a priori estimate (2.42), the expected rate of convergence is po + | if 
the solution is smooth in K € Vo. To prohibit p-enrichments in discontinuous 
regions, we simply set p L = p 0 ,L = 1 ,n K . if fi K < p 0 + | for the parent 
element. The contribution of these elements to the global error therefore re- 
mains fixed in the p-step of the adaptive strategy. If this contribution exceeds 
the target Ox, then an additional /i-step is required before the p-step in order 
to achieve the target error. 


, K = l,-.-,N(V 0 ) (4.15) 
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4.3 An hp - step as an Alternative to the p - Step 

For smooth solutions or for solutions which contain mild discontinuities, the 
p-step should be adequate to reduce the global error to the specified level. For 
problems with strong discontinuities, however, the local error at discontinuities 
may be significantly large and dominate the global error obtained after the h- 
step, particularly in the current implementation where a maximum of 2 levels 
of refinement are permitted. 

An alternative to the p-step of the adaptive strategy is an hp- step where 
/i-refinement is performed at the discontinuity and p-enrichment is performed 
in smooth regions. In this case, the target error is specified as a reduction factor 
for the error in discontinuous regions and for the smooth region individually. 
Let r] D denote the normalized error in discontinuous regions and tj s denote the 
normalized error in smooth regions. Then the target error for the hp- step is 
specified by the reduction factors a D and a s as 

V t = \/( a D 7 ?r,) 2 + (^V/s) 2 ( 4 - 16 ) 

The criteria of p K < p 0 + |, where p K is given by (4.15), is used to distin- 
guish discontinuous regions from smooth regions. Equations (4.8) and (4.9), 
with N(Vq) replaced by the number of elements in the discontinuous region, 
are used to determine the number of elements required to reduce the error in 
discontinuous regions to a D r] D . Equation (4.14), with N(V i) replaced by the 
number of elements in the smooth region, is used to determine the values of 
p K required to reduce the error in the smooth region to a s Tj s . This approach 
is referred to as an hhp- adaptive strategy. 



53 


4.4 Selection of the Parameters 

Formally, the parameters p K and v K depend on the global regularity of the 
solution. However, the rate of convergence of the local error depends on the 
local regularity of the solution. For piecewise continuous solutions, the rate of 
convergence of the local error varies greatly between the smooth regions of the 
solution and some small neighborhood around discontinuities. Using global val- 
ues of these parameters based on an irregular solution results in over-refinement 
of smooth regions while using global values based on smooth solutions results 
in under-refinement of discontinuities. Here we use local values of p K and 
v k which are initially computed for a uniform /i-refinement and a uniform p- 
enrichment of a coarse mesh. These local values are passed onto the initial 
mesh used for the adaptive strategy. Local values of p K are then re-computed 
after the h - step using (4.15). These are the values used in (4.12) for the p- 
or /ip-step. While there is little theoretical justification for using local values, 
numerical results indicate that the approach works quite well for solutions with 
discontinuities. 

Selection of the reduction factor a used to determine the intermediate 
target error for the A-step of the adaptive strategy is important in obtaining an 
optimal mesh. Specifying a value of a which gives an intermediate error which 
is closer to the target error than it is to the initial error will result in meshes 
with mostly ^-refinement. Specifying a value of a which gives an intermediate 
target error which is closer to the initial error than it is to the target error 
leads to meshes with little /i-refinement and elements with large values of p K . 
Numerical experiments for elliptic problems suggest that the optimal choice is 
a = j;, where p and v are the rates of convergence of the global error. [41] 


Chapter 5 


Numerical Examples 


The discontinuous Galerkin method is used to solve several examples to verify 
the a priori error estimates derived in Chapter 2, and to investigate the perfor- 
mance of the the a posteriori error estimates of Chapter 3 and the ftp-adaptive 
strategy of Chapter 4. The reliability of the a posteriori error estimates is mea- 
sured by the effectivity index which is the ratio of the estimated error to the 
exact error. A reliable estimate is one for which the efFectivity index is close to 
one. 


5.1 Example 1 

We solve the linear model problem (2.1) with the following data: 

(i) n = C-i,i) x (—i,i) 

(ii) (3 = (0.8, 0.6) r 

(iii) a(x) = 1.0 

(iv) g = 1.0 
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Figure 5.1: Quadrilateral element mesh used for quasi-uniform refinements. 

The source term / is chosen so that the exact solution to (2.1) is the 
function, 

u(x,y) = 1 + sin (|(1 + x)(l + J/) 2 ) (5.1) 

The a priori error estimate (2.42) is verified by solving the problem for a se- 
quence of uniform /i-refinements and p-enrichments of a mesh of square elements 
and quasi-uniform /i-refinements and p-enrichments of the mesh of quadrilat- 
eral elements shown in Fig. 5.1. The mesh-dependent norm of the actual error 
in the solution obtained with varying h and p is listed in Table 5.1 for the 
square element mesh and in Table 5.2 for the quadrilateral element mesh. 
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Mesh 

— log h 

log |||t/ Ufa 1 1 \hp t l3 

p = 1 

p = 2 

P = 3 

p = 4 

2x2 

0.000 



1.8323 

2.2787 

4x4 

0.301 

0.5552 

1.7066 

2.5426 

3.6065 

8x8 

0.602 

0.9692 

2.3909 

3.5467 

4.9612 

16 x 16 

0.903 

1.4003 

3.1163 

4.5834 

6.3047 

32 x 32 

1.204 

1.8412 

3.8574 




Table 5.1: Example 1 - Error using uniform hp meshes. 


Mesh 

log h 

1 

£ 

1 

e 
1 ;*■ 

p = i 

P = 2 

P = 3 

p = 4 

2x2 

-0.2116 


0.8586 

1.7402 

2.2831 

4x4 

0.0689 

0.5153 

1.5930 

2.5395 

3.4998 

8x8 

0.347 

0.9571 

2.3641 

3.5814 

4.9723 

16 x 16 

0.641 

1.3913 

3.0955 

4.6208 

6.3196 

32 x 32 

0.938 

1.8129 

3.7870 

i 



Table 5.2: Example 1 - Error using quasi-uniform h and uniform p meshes. 
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Figure 5.2: Example 1- Rate of convergence of error for fixed p. 

To verify the estimate (2.42), first consider the case when p K is fixed 
and h K is varied. According to (2.42), |||e||| Api/ j < Ch P K K+ * ||«|| r>n . This is 
verified in Fig. 5.2 where |||e|||^ p ,/j is shown as a function of h K . On the log- 
log scale, the slope of the lines corresponding to a fixed value of p K is indeed 
P K + I f° r both the uniform and quasi-uniform meshes. Next consider the case 
when h K is fixed and p K is varied. In this case, the estimate (2.42) reduces to 
lll e IIUp,/0 ^ C f p Jf _r+1 |M| r) n. Since u € C°°(f2), exponential rates of convergence 
are expected. This is confirmed in Fig. 5.3 where the curves corresponding to 
lll e IIUp,/3 as a function of p K have a slope on the log-log scale which increases 
as p K increases. These results are combined in Fig. 5.4 where |||c|||fc PiJ 9 is 
shown as a function of the total number of unknowns in the solution. The 
solid lines represent /i-refinements for a fixed p and the dashed lines represent 
p-enrichment for a fixed h. Clearly for smooth solutions, higher-order accuracy 
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Figure 5.3: Example 1- Rate of convergence of error for fixed h. 
is achieved for the same number of unknowns using higher-order elements. 

Next we investigate the performance of the a posteriori error estimates 
in Chapter 3. Recall that the complete polynomial space Q Pk+ ° k (/if) or the 
incomplete space Q h A (/if) \ Q K (/if), a K > 1 may be used in approximating 
the solution to the local problem. The effect of approximating the local problem 
on the performance of the error estimate for the lower bound (3.1) is shown in 
Table 5.3. The effectivity indices listed in Table 5.3 are greater than one for all 
values of a K when the complete polynomial space is used and less than one for 
small values of a K when the incomplete polynomial space is used. Note that 
the effectivity index is closest to one for the complete polynomial space with 
cr K = 1 . The effect of approximating the local problem on the performance 
of the error estimate for the upper bound (3.7) is shown in Table 5.4. The 
effectivity indices for the upper bound estimate are significantly larger than 



<Pk € Q Pk+<7 '<(!<) 

<P k €Q^k(K)\QPk( K ) 

Mesh 

P JL 

a K 

Vi 

Vl 

8x8 

1 

1 

1.0938 

0.7628 

8x8 

1 

2 

1.1272 

0.8921 

8x8 

1 

3 

1.1347 

0.9788 

8x8 

1 

4 

1.1372 

1.0175 

16 x 16 

1 

1 

1.1765 

0.7933 

16 x 16 

1 

2 

1.2229 

0.9492 

16 x 16 

1 

3 

1.2340 

1.0465 

16 x 16 

1 

4 

1.2378 

1.0898 

8x8 

2 

1 

1.1224 

0.9555 

8x8 

2 

2 

1.2009 

1.0711 


Table 5.3: Example 1 - Effect of the approximation of the local problem for 
the lower bound on the effectivity index. 
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Figure 5.4: Example 1- Rate of convergence of error with number of unknowns. 

one when the complete polynomial space is used and close to one when the 
incomplete polynomial space with a K = 1 is used to approximate the solution 
to the local problem. 

Next we verify that the error estimate exhibits the same rates of con- 
vergence as the actual error with /i-refmement or p-enrichment. Based on 
the results from above, the error is estimated by solving the lower bound lo- 
cal problem in Q (A) and by solving the upper bound local problem in 
Q h (K) \ Q h (K). The estimated error for a sequence of ^-refinements of 
meshes with fixed p is shown as a function of the mesh size in Fig. 5.5. The 
slope of the lines is p + | as in the case of the actual error (see Fig. 5.2). The 
estimated error for a sequence of p-enrichments of a uniform mesh is shown in 
Fig. 5.6 where the same behavior as the actual error (see Fig. 5.3) is observed. 




*1> K € Q p «+°1< ( K ) $ K e Q p >< + °‘< (I<) \ Q p ‘< (K) 


Mesh 

Pk 


Pu 

Vu 

8x8 

1 

1 

4.1911 

1.0536 

8x8 

1 

2 

4.3603 

1.3551 

8x8 

1 

3 

5.1297 

2.0930 


8x8 

1 

4 

5.3218 

2.6803 


16 x 16 

1 

1 

6.4441 

1.1238 


16 x 16 

1 

2 

6.6729 

1.4980 


16 x 16 

1 

3 

7.9048 

2.8298 


16 x 16 

1 

4 

8.1157 

3.6008 


8x8 

2 

1 

2.6353 

1.1957 


8x8 

2 

2 

3.7052 

1.4012 



Table 5.4: Example 1 - Effect of the approximation of the local problem for 
the upper bound on the effectivity index. 






64 


1.07 

1.02 

0.97 

0.92 

0.88 

0.83 

0.78 

0.73 


Figure 5.7: Example 1 - Local effectivity index for error estimate based on the 
upper bound local problem. (8x8 mesh, p = 1) 

While the theory developed thusfar applies to global error estimates, 
local effectivity indices near unity are desired in order to use the estimate to 
drive an effective adaptive strategy. The local (element) effectivity index for 
the error estimate based on the upper bound local problem (3.7) using the 
incomplete polynomial space Q h (K) \ Q K (K) is shown for a uniform 8x8 
element p = 1 mesh in Fig. 5.7, for a uniform 8x8 element p — 2 mesh in Fig. 
5.8, and for a uniform 16 x 16 element p — 1 mesh in Fig. 5.9. For all cases 
investigated, the local effectivity index is close to one except in a few isolated 
elements. This indicates that the error estimate is reliable enough to drive the 
/ip-adaptive strategy. Therefore, the solution to the upper bound local problem 
on the incomplete polynomial space with a K — 1 will be used throughout to 
drive the /tp-adaptive strategy. 
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1.25 

1.22 

1.19 

1.17 

1.14 

1.11 

1.08 

1.05 


Figure 5.8: Example 1 - Local effectivity index for error estimate based on the 
upper bound local problem. (8x8 mesh, p = 2) 



1.11 

1.04 

0.97 

0.90 

0.84 

0.77 

0.70 

0.64 


Figure 5.9: Example 1 - Local effectivity index for error estimate based on the 
upper bound local problem. (16 x 16 mesh, p = 1) 
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Estimated Error 



9.3E-2 
8.2E-2 
7.1 E-2 
6.0E-2 
4.9E-2 
3.8E-2 
2.8E-2 
1.7E-2 


Figure 5.10: Example 1 - Error distribution on initial mesh. 

The results of applying the adaptive strategy described in Chapter 4 
to this problem are summarized in Table 5.5. The normalized error listed in 
the table is the ratio of the global error to the sum of the global solution and 
the error in the norm associated with the local problem. Starting with the 
estimated error on an initial mesh of 4 x 4 p = 1 elements, shown in Fig. 5.10, 
a target error of 1.5 percent is specified. An estimated error of 1.2 percent is 
actually achieved on the resulting h - adapted mesh (see Fig. 5.11). For the 
p-step of the adaptive strategy, a target error of 0.1 percent is specified. An 
estimated error of 0.11 percent is actually achieved on the p - adapted mesh (see 
Fig. 5.12). The distribution of the error on the adapted meshes, shown in Figs. 
5.11 and 5.12, is also reduced by an order of magnitude at each adaptive step. 


5 

4 

3 

2 

1 


1.4E-2 

1.2E-2 

1.0E-2 

8.7E-3 

7.0E-3 

5.3E-3 

3.6E-3 

1.9E-3 
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8.8E-4 

7.7E-4 

6.6E-4 

5.5E-4 

4.4E-4 

3.3E-4 

2.2E-4 

1.1E-4 


Figure 5.12: Example 1 - Mesh and error distribution after the p-step. 


5.2 Example 2 

We solve the linear model problem (2.1) with the following data: 


69 


(i) ft = (-1,1) x (-1,1) 

(ii) {3 = (1.0,0.0) T 


(iii) a(x) = 1.0 

r \ _ / 3e~ HlW) if y < 0 

? \ —3e -5 ( 1+y2 ) otherwise 


The source term / is chosen so that the exact solution to (2.1) is the discon- 
tinuous function, 


f 3e- 5 (* 2 +jd) if 3/ < 0 
( — 3e _5 ^ 2+y2 l otherwise 


(5.2) 


The discontinuity is aligned with element interfaces at y = 0 to illustrate the 
advantage of using a discontinuous method to capture discontinuities, particu- 
larly if the adaptive scheme includes some shock fitting which aligns the grid 
with the discontinuity. 


The problem was solved using a variety of uniform meshes with ft- 
refinements, p-enrichments, and the /ip-adaptive strategy with no special treat- 
ment at the shock. The error histories for two ftp-adaptive solutions with dif- 
ferent initial meshes are listed in Tables 5.6 and 5.7. For both cases, the target 
error was achieved at each step in the adaptive process. Note also that the 
global effectivity index is near unity for all steps. 

The rate of convergence of the estimated and exact error is compared 
in Fig. 5.13. The exact error (denoted by a solid line in the figure) and the 
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Adaptive 

step 

Target 

normalized 

error 

Exact 

error 

IMLp 

Estimated 

error 

m Au 

Achieved 

normalized 

estimated 

error 
I M\a„ 

Initial 4 x 4 p — 1 mesh 


0.1663 

0.2078 

0.03786 

h-refinement 

0.015 

0.03619 

0.03847 

0.012105 

p-enrichment 

0.001 

0.00377 

0.00355 

0.001116 


Table 5.5: Example 1 - Error history for an adaptive hp solution. 


Adaptive 

step 

Target 

normalized 

error 

Achieved 

normalized 

error 

Estimated 

error 

\m\ Au 

Effectivity 

index 

Initial 4x4 mesh, p = 2 


0.091 

0.371 

1.055 

/i-refinement 

0.05 

0.031 

0.127 

1.073 

p-enrichment 

0.005 

0.0029 

0.012 

1.424 


Table 5.6: Example 2 - Error history for an adaptive hp solution starting from 
a uniform 4x4 mesh, p = 2. 
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Adaptive 

step 

Target 

normalized 

error 

Achieved 

normalized 

error 

Estimated 

error 

M\a 0 

Effectivity 

index 

Initial 8x8 mesh, p = 1 


0.154 

0.616 

0.998 

ft-refinement 

0.075 

0.033 

0.137 

0.996 

p-enrichment 

0.005 

0.0055 

0.023 

0.901 


Table 5.7: Example 2 - Error history for an adaptive ftp solution starting with 
a uniform 8x8 mesh, p = 1. 


estimated error (denoted by a dashed line) are in close agreement, indicating 
the reliability of the estimate. Note that with the discontinuity aligned with 
element interfaces, the error behaves as if the solution is smooth, that is, alge- 
braic rates of convergence are achieved with respect to mesh refinements, and 
exponential rates of convergence are achieved with respect to p-enrichments. 
In this case, the most significant error reduction with fewest degrees of freedom 
will result by specifying a target error for the ft- step which is closer to the initial 
error than to the final target error. This is verified by the two curves corre- 
sponding to the ftp-adaptive solutions in Fig. 5.13. Results of the ftp- adaptive 
solution from the initial 8x8 element p = 1 mesh are shown in Figs. 5.14 - 
5.20. The estimated error in the solution on the initial mesh is shown in Fig. 
5.14 with the corresponding effectivity index shown in Fig. 5.15. 

For the ft-step in the adaptive procedure, a normalized target error of 
7.5 percent resulted in the mesh shown in Fig. 5.16. The estimated error in 
the solution obtained on the ft-adapted mesh is also shown in Fig. 5.16 and 
the corresponding effectivity index is shown in 5.17. Poor local error estimates 
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Figure 5.13: Example 2 - Rate of convergence of the error with respect to the 
total number of unknowns. 

are observed in the two parallel vertical regions indicated by the darker shades 
in Fig. 5.17. Moreover, the local error is significantly underestimated in these 
regions, possibly due to a failure of the procedure to adequately handle the 
very high changes in gradients in these regions. The global efFectivity indices, 
however, are quite satisfactory with efFectivity indices very near unity. 

For the p-step, a normalized target error of 0.5 percent resulted in the 
distribution of p shown in Fig. 5.18. The estimated error for the solution 
obtained on the ftp-mesh and corresponding local efFectivity index are shown 
in Figs. 5.19 and 5.20. The same degradation of the local efFectivity indices are 
observed in the regions where up is large. Moreover, there is a slight decrease in 
the global efFectivity index, however, the global value is still quite acceptable. 
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0.0257 

0.0221 

0.0185 

0.0149 

0.0113 

0.0077 

0.0041 

0.0005 


Figure 5.16: Example 2 - Estimated error on ft-adapted mesh. 



1.0600 

0.9414 

0.8229 

0.7043 

0.5857 

0.4671 

0.3486 

0.2300 


Figure 5.17: Example 2 - Local effectivity index for error estimate on /i-adapted 
mesh. 



75 



In 4.0000 

| :ari 3.0000 
2.0000 
1 .0000 


Figure 5.18: Example 2 - Adaptive p-enriched mesh. 



0.0035 

0.0031 

0.0026 

0.0022 

0.0018 

0.0013 

0.0009 

0.0004 


Figure 5.19: Example 2 - Estimated error on adaptive p-enriched mesh. 
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1.6377 
1.4361 
1 .2346 
1.0330 
0.8315 
0.6299 
0.4284 
0.2269 


Figure 5 . 20 : Example 2 - Local effectivity index on adaptive p-enriched mesh. 


5.3 Example 3 

The following data is used in (2.1): 

(i) ft = (-1,1) x (-1,1) 

(ii) ,3 = (f , f ) T 

(iii) a(x) = 1.0 


(iv) g(x,y ) = 


5 e U +v 1 + 3 e I 1 "*'! 1 ' 2FI x = — 1 
-1 - 8e- 5 K x ~i> 2+ i] y = - 1 


The source term / in (2.1) is chosen so that the exact solution is a function 
which is discontinuous along the domain diagonal given by 


f 5 e-K*H) 2 +y 2 ] + 3 e-(^+( 3 /-?) 2 l if y > x 
\ — 1 — 8e -5 K r ~2) 2+ ( !/+ 2) 2 l otherwise 
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The global effectivity index for the estimate obtained by solving the 
upper bound local problem in the space Q p k +1 (K)\Q p k(K ) for several uniform 
hp meshes is listed in Table 5.8 which shows that the global error is slightly 
under-estimated. Results at each step in the adaptive strategy, shown in Figs. 
5.21 - 5.27, show that the under-estimation of the global error is primarily due 
to the under-estimation of the local error at the discontinuity. A summary of 
the hp-strategy is listed in Table 5.9. Note that while the adaptive strategy is 
able to reduce the error to the target value in the h-step, the achieved error 
after the p-step largely represents the remaining error in the discontinuity after 
the h-step. 

The error achieved by the adaptive strategy is compared to uniform 
refinements of p = 1 and p = 2 meshes in Fig. 5.28. This plot shows that 
the rate at which the error is reduced in the p-step is much higher than is 
possible using h — refinements alone. Moreover, the main source of the error in 
the final solution is attributed to the discontinuity (see Fig. 5.26. In Chapter 
6, we solve a similar problem on a parallel computer and use the hhp - adaptive 
strategy described in Chapter 4 to reduce the error in solution over the entire 
domain. 


5.4 Example 4 

The following data is used for (2.1): 

(i) fl is an equilateral triangle with side length of 2 and the base al- 
ligned with the x-axis 


(ii)/3 = (f,f ) r 
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Mesh 

P,< 

1 \M\av 

IMU* 

% 

4x4 

1 

3.01325 

3.81716 

0.79 

4x4 

2 

1.49536 

1.97645 

0.76 

8x8 

1 

1.62507 

2.10238 

0.77 

8x8 

2 

0.82209 

1.35611 

0.61 

16 x 16 

1 

0.88158 

1.37317 

0.64 

16 x 16 

2 

0.63044 

1.00438 

0.63 

32 x 32 

1 

0.85472 

1.01102 

0.85 


Table 5.8: Example 3 - Error estimate obtained by approximating the upper 
bound local problem in Q p k + 1 (I\) \ Q p i<(K). 


Adaptive 

step 

Target 

normalized 

error 

Achieved 

normalized 

error 

Estimated 

error 

m Au 

Effectivity 

index 

Initial 8x8 element p = 1 mesh 

— 

0.118 

1.625 

0.77 

h - refinement 

0.05 

0.048 

0.673 

0.63 

p-enrichment 

0.03 

0.038 

0.541 

0.55 


Table 5.9: Example 3 - Error history for an adaptive hp solution starting from 
a uniform 8x8 mesh, p = 1. 














0.1032 

0.0903 

0.0774 

0.0645 

0.0516 

0.0388 

0.0259 

0.0130 


Figure 5.23: Example 3 - Estimated error on /i-adapted mesh. 




1.67 

1.48 

1.29 

1.11 

0.92 

0.73 

0.54 

0.35 


Figure 5.24: Example 3 - Local effectivity index on /i-adapted mesh. 
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(iii) a(x) = 1.0 

The source term / in (2.1) is chosen so that the exact solution is a 

u(r) = sinSrjtan* 1 100 + tan _1 [100(r — 1)]} (5.4) 

where 

r 2 = (x + 0.1) 2 + Ay 2 

The inflow boundary conditions are obtained by evaluating the exact solution 
along the inflow boundary 

, , __ ( u(r) r 2 = (x -f 0T) 2 on y = 0 

| w(r) r 2 = (x + 0.1) 2 + 4 tan 2 (|)x 2 on y = tan(|)x 

The error history for an hp- adaptive solution is listed in Table 5.10 and the 
resulting meshes, local error estimates, and local efFectivity indices are shown 
in Figs. 5.29 - 5.35. 

While the exact solution to this problem is continuous, it contains a 
very steep front which can be seen as the dark regions in Fig. 5.29. The global 
efFectivity indices listed in Table 5.10 demonstrate the reliablility of the error 
estimate on non-rectangular meshes. The local efFectivity indices for all the 
meshes are close to unity over most elements, however, there is some slight 
under-estimation of the error for the uniform p-meshes. Though it is difficult 
to see from the figure, the under-estimation of the error occurs in the region of 
the steep front. There is also some rather severe over-estimation of the error 
on the p-adapted mesh. This over-estimation occurs primarily in the p = 5 
elements and does not have a significant effect on the global efFectivity of the 
error estimate. 
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Adaptive 

step 

Target 

normalized 

error 

Achieved 

normalized 

error 

Estimated 

error 

uk 

Effectivity 

index 

Initial mesh p = 1 


0.130 

2.514 

0.87 

^-refinement 

0.075 

0.029 

0.641 

1.027 

p-enrichment 

0.003 

0.0025 

0.559 

1.028 


Table 5.10: Example 4 - Error history for an adaptive hp solution. 



0.6833 

0.5979 

0.5125 

0.4271 

0.3417 

0.2562 

0.1708 

0.0854 


Figure 5.29: Example 4 - Estimated error on initial mesh, p = 1. 
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1.36 

1.21 

1.07 

0.92 

0.77 

0.62 

0.48 

0.33 


Figure 5.30: Example 4 - Local effectivity index on initial mesh, p = 1. 



0.0960 

0.0840 

0.0720 

0.0600 

0.0480 

0.0360 

0.0240 

0.0120 


Figure 5.31: Example 4 - Estimated error on /i-adapted mesh 
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1.36 

1.21 

1.07 

0.92 

0.77 

0.62 

0.48 

0.33 


Figure 5.32: Example 4 -Local effectivity index on h - adapted mesh. 



MmMw. 

iNttlk 



5.00 

4.00 

3.00 

2.00 
1.00 


Figure 5.33: Example 4 - p-adapted mesh. 




0.1027 

0.0898 

0.0770 

0.0642 

0.0513 

0.0385 

0.0257 

0.0128 


Figure 5.34: Example 4 - Estimated error on p-adapted mesh. 



10.12 

8.86 

7.61 

6.35 

5.09 

3.84 

2.58 

1.33 


Figure 5.35: Example 4 - Local effectivity index on p-adapted mesh. 


Chapter 6 


Parallel Implementation 


The time-marching versions of the discontinuous Galerkin methods, (2.54), 
(2.55), and the Runge-Kutta discontinuous Galerkin method described in the 
next Chapter, fall naturally into the class of single program multiple data 
(SPMD) parallel applications. Given the solution at time level t n , the solution 
is advanced to time level < n+i by solving a ( p K + l) 2 x ( p K + l) 2 system of 
linear equations for the ( p K + l) 2 degrees of freedom for every element K in the 
partition. The only coupling between elements in the partition arises in the 
evaluation of the boundary integrals in (2.8) where the solution along common 
edges of neighboring elements is needed. The evolution of the solution can be 
performed on all the elements simultaneously, once this information is available. 

The primary issue in a parallel implementation of discontinuous Galerkin 
methods is to balance the workload among the available processors while min- 
imizing the communication between processors, thereby optimizing the utiliza- 
tion of the multi-processor environment. For a machine with P processors, 
this is accomplished by dividing the partition into P subdomains and assigning 
the elements contained in a particular subdomain to a particular processor. 
In order to minimize communications, the interface of the subdomain bound- 
aries should have as small a measure as possible. Moreover, the local nature 
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of the discontinuous Galerkin formulation requires communication of the solu- 
tion only along subdomain boundaries at each time step. Since computations 
are performed at an element level, communications can be overlapped with 
the computation to further minimize the penalty of communicating between 
processors. 

While the parallel implementation described in this Chapter is targeted 
for the Intel iPSC 860 computer which is a distributed memory machine with 
32 processors arranged in a hypercube architecture, many concepts are general, 
and, therefore applicable to other multi-processor machines. 

6.1 Domain Decomposition for hp Meshes 

The goal of the domain decomposition strategy is to evenly distribute the work- 
load among the processors while minimizing the size of the subdomain bound- 
aries. Most domain decomposition methods have been developed and analyzed 
for h - type meshes where the number of degrees of freedom, and hence the com- 
putational effort, is the same for every element in the mesh. For these types of 
meshes, equally distributing the elements among the available processors will 
result in a balanced load. 

The most successful domain decomposition methods in this situation 
are based on recursive bisectioning of either the physical domain or an order- 
ing of the elements. In the recursive bisectioning of the physical domain, trial 
separators define possible subdomain configurations. The selection of a sepa- 
rator as a subdomain interface is based on the resulting load balance as well as 
interface size. Vavasis [52] has obtained theoretical bounds on the achievable 
load balance and interface size. One disadvantage of this approach is that it 


< 2 .- 2 ,. 
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can be computationally expensive for multiple space dimensions. 

Recursive bisectioning methods based on an ordering of the elements 
have a computational advantage since the bisectioning is performed on a one- 
dimensional list of elements, regardless of the spatial dimension of the domain. 
One difficulty with this approach is constructing an ordering which preserves 
the locality of the elements in the mesh. A locality-preserving ordering is neces- 
sary to avoid multiply connected or disconnected subdomains and to minimize 
interface size. Pot.hen, Simon, and Liou [45] construct such an ordering by 
using the second eigenvector of the Laplacian matrix associated with the graph 
of the mesh. 

For hp meshes, where the number of degrees of freedom (and hence the 
computational effort) vary from element to element, the domain decomposition 
must include some measure of the computational work for each element. Here 
we investigate two load-based recursive bisection methods developed by Patra 
[40] which seek to balance the load using some local measure of the computa- 
tional work. In this study, the number of degrees of freedom in an element, 
(p + l) 2 , are used as a measure of its computational load. 

6.1.1 Recursive Load Based Bisection of Coordinates (RLBBC) 

An algorithm for this method which is based on recursive bisectioning of the 
physical domain is given below. 

6k'- computational effort estimate for element K, 6k specified as dof in 
the description of the algorithm. It may be replaced by any alternate measure 
of computational effort. 
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Dj : list of elements in subdomain I. 
ni : number of trial separator surfaces. 
qi : quality index for a trial separator. 

1. Compute maximum and minimum coordinates in any one of the dimensions 

of the entire domain rl m ^ n , xl max 

For i = 1 to n/ do 

2. compute 

™ | . rp 1 I (^moi-Uffun) 

sL kL ± mtn T 

n x 

( l l ~ do fright * do f tot 4 " dofi nter 

where dof le j t and do f right are the degrees of freedom to the left and right 
of xl i, respectively, and dof inter is the degrees of freedom on the trial 
separator xl,. 

3. Choose as interface the separator corresponding to the lowest qj. 

4. If the center of mass of element has an xl coordinate less than that of the 

interface then 

Dy 4- D x U {I<} 

Else 

D 2 «- D 2 U {K} 

At this stage, the original domain has been split into two. 

5. For the next level of decomposition apply steps 1-4 with Dy and D 2 instead 

of the entire domain and with x 2 as the coordinate. This will result in 
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four subdomains. In three-dimensions, for the next level use 1-4 with 
these four subdomains and with X 3 as the coordinate. This process is 
recursively continued until the desired number of domains is attained. 
Clearly, for better shaped domains, equal numbers of splits in each coor- 
dinate must be made. Thus, for two dimensions, 4 n subdomains and in 
three dimensions, 8 n subdomains are obtained. 

6.1.2 Recursive Load Based Bisection of an Ordering (RLBBO) 

This method is based on recursive bisectioning of an ordering of the elements. 
The ordering is based on results of Peano[42] and Hilbert [27] concerning a 
class of continuous mappings of the unit interval onto a unit hypercube. The 
significance of their results is that one can construct a space filling curve which 
connects a set of points in rc-dimensional space and uniqely maps them onto the 
unit interval. Applying this mapping to the set of points given by the element 
centroids, results in an ordering of the elements defined by its location on the 
unit interval. Complete details of the Peano-Hilert ordering and proof that it 
is locality-preserving can be found in Patra [40]. 

An algoritm for this method follows: 

h n : 3ft — ► U n , where Ui is a the unit interval and U n is the unit hyper- 

cube. 

D{ : list of elements in subdomain I. 

rii : number of trial separator surfaces. 

qi : quality index for a trial separator. 


1. Create an ordering of the elements by mapping the centroids of the ele- 
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ments onto a Peano-Hilbert curve. (See [40]) 

2. Let tx be the distance of the centroid of element K along the space filling 
curve. 

3. Compute maximum and minimum of tx, t m ax and t m in. 

4. Compute n/ trial separator levels as 

^ 4 i tmax tmin 

t't — t'min i 

Til 


5. For each t t compute a quality of interface index qi 

q t = abs(— - - 1) • do f tot -f dofi nter 

fright 

Replace dof by error or other load estimate as appropriate. 

6. Choose as interface t{ nt the t{ that corresponds to lowest qi 


7. For K = 1 to the total number of elements 


If t K < ti nt then 
D x f-AU {I<} 
else 

D 2 * — D 2 U { A } 
end if 
end for 

At this stage, the original domain has been split into two. 

8. Apply 1-7 recursively on each of the generated subdomains. 
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6.2 Communications 

Communication between processors on distributed memory machines can sig- 
nificantly effect the overall performance of the parallel implementation, partic- 
ularly if a processor must wait to receive information from another processor 
before proceeding with a calculation. To compute the solution at time level t n +\ 
on a particular subdomain, the solution at t n is needed from interface elements 
on neighboring subdomains. Using synchronous communications, that is, re- 
questing information from neighboring processors at the time that it is needed, 
leads to unnecessary waiting by all processors. Moreover, communication con- 
flicts are likely to occur since two-way communication is required across interior 
subdomain boundaries. Asynchronous communications are used to minimize 
this wait time. When the solution for an interface element is computed on a 
processor, it is then sent to the processor containing the neighboring element. 


6.3 Numerical Results 


Results are presented for an Mp-adaptive solution to an example similar to 
example 3 in Chapter 5, except that the source term / is selected so that the 
exact solution is 


f 5 e -E(^+|) 2 +y 2 ] + 3 e -U 2 +(i^) 2 l if y > x 
\ —1 — 5 e~ 5 E( I- j) 2 +(v+ 2 ) 2 ] otherwise 


( 6 . 1 ) 


The error history obtained with the hhp - adaptive strategy described in 
Chapter 4 is listed in Table 6.1 which shows that this approach is very effective 
at reducing the global error in the discontinuous solution. These results are 
compared graphically with results from Chapter 3 in Fig. 6.1. Note the sig- 
nificant decrease in the error due to the additional ^-refinement at the shock. 
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Adaptive 

step 

Target 

normalized 

error 

Achieved 

normalized 

error 

Estimated 

error 

IWL, 

Effectivity 

index 

Initial 8x8 mesh, p = 1 


0.126 

1.391 

0.72 

/^-refinement 

0.0628 

0.0759 

0.861 

0.69 

p-enrichment 

0.0558 

0.0311 

0.359 

0.53 


Table 6.1: Error history for an ftftp- adaptive solution starting from a uniform 
8x8 mesh, p = 1. 


While this adaptive strategy leads to significantly more unknowns in the prob- 
lem, the reduction in the error, when compared to error obtained by uniform 
ft-refinements of a mesh of p = 1 or p = 2 elements, is orders of magnitude. 
The final hp mesh, shown in Fig. 6.2, is particularly demanding for a domain 
decomposition method because of the concentration of elements resulting from 
4 levels of refinement in the lower left quadrant of the domain and the large p 
elements in the smooth regions. 

We use the speedup as a measure of the effectiveness of the domain 
decomposition startegy at balancing the load. The speedup is defined here as 
the ratio of the CPU time required to obtain the solution using 4 processors 
to the CPU time required to obtain the solution using n processors, where 
n > 4. Ideal (also called linear) speedup implies that the work load is equally 
balanced among the processors and is indicated by a slope of two on a graph of 
the speedup as a function of the number of processors. The speedup obtained 
on the ftp-mesh using the two decomposition strategies is shown in Fig. 6.3. 
Nearly optimal speedup of 1.8 is obtained using the RLBBC method when going 
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number of unknowns 

Figure 6.1: Rate of convergence of the error with respect to the total number 
of unknowns. 

from 4 to 8 processors. The speedup drops off as the number of processors 
increases, indicating that the load is unbalanced. The speedup of 1.3 when 
going from 4 to 8 processors when using the RLBBO method indicates that 
the method yields very poor load balancing when only a few processors are 
used. Fortunately, the RLBBO method provides more balanced loads for a 
larger number of processors where a speedup of 1.85 is achieved when going 
from 8 to 16 processors. 









Chapter 7 


Extensions to Nonlinear Hyperbolic 
Conservation Laws 


In this chapter, we describe extensions of the discontinuous Galerkin method to 
nonlinear hyperbolic conservation laws of the form (1.5) on convex polygonal 
domains. For simplicity, the method is described for a two-dimensional scalar 
conservation law written as 

u t + V • q(u) = 0 (x, y) e ft C 3ft 2 , t > 0 (7.1) 

where q (u) — f l {u) i-f f 2 (u) j denotes the flux vector. For a description of the 
method for hyperbolic systems of conservation laws, in particular, the Euler 
equations of gas dynamics, see Bey and 0den[6]. 

To complete the initial boundary- value problem, (7.1) is accompanied 
by initial conditions 

w 0*h2/,0) = u 0 (x,y); (x,y) e fl (7.2) 

and boundary conditions of the form 

u(x,y,t) = d{t)] (x,y) (7.3) 

where T - is the inflow boundary, F - = {(a;,?/) £ dfi : • n < 0} , n is the 

outward unit normal to the boundary (dfl = T - fl T + ), and d(t) is prescribed 
inflow data. 
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The version of the discontinuous Galerkin method used here, first pro- 
posed by Cockburn and Shu [13], [12], is based on the method of lines. The sys- 
tem of ordinary differential equations resulting from the discontinuous Galerkin 
spatial approximation is marched in time using a multi-stage Runge-Kutta 
scheme. A local projection is used at each stage to control oscillations and 
prevent nonlinear instabilities. In reference [13], a projection is constructed 
for p — 1 and p — 2 elements in one space dimension and possible extensions 
to higher-order elements are proposed. In reference [12], some important the- 
oretical results are derived for multi-dimension problems, but a projection is 
constructed and verified only for linear p = 1 triangles. Moreover, the emphasis 
of their work is in predicting element mean values, and while the high-order 
solution is available, it is ignored in the presentation of numerical results. Here 
we construct a very simple projection for general quadrilateral (or triangular) 
elements of arbitrary order and take full advantage of the finite element approx- 
imation. For alternate projection strategies for quadrilaterals, we refer to Bey 
and Oden [6] and to Biswas, Devine, and Flaherty [7] for a general projection 
strategy on Cartesian grids. 

7.1 Discontinuous Galerkin Spatial Approximation 

The primary difference in the discontinuous Galerkin formulation for nonlinear 
conservation laws is the treatment of the boundary fluxes at element interfaces 
and the need for the projection to control oscillations. Recall that for the linear 
conservation law q (u) = flu and hence ^ = fl • < 0, which defines 

element inflow boundaries, is known a priori. For the nonlinear conservation 
law, ^ • riA' is a function of u and is thus unknown. 
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The discontinuous Galerkin spatial approximation is defined on a par- 
tition of the domain satisfying the properties in Chapter 2 (except property 
(v) defining element inflow boundaries) and solutions are sought in the finite 
dimensional space V p (Vh ) = {v £ £ a (f2) | i>| £ Q p x(K)}. The weak form of 

(7.1) on a single element is obtained by multiplying by a test function u(x) and 
integrating the first-order spatial term by parts: 


Find u(x, J)| A . £ Zd(/v) x C l (0, T) such that 


j u(x, 0)v(x)dx = J u 0 (x)v(x)dx and 
— (u, v) — I q(u) • Vvdx + f q(u) • n kv ds — 0 

ut J I\ J 3 hi 

for all admissible test functions i?(x). 


(7.4) 

(7.5) 


The element boundary flux q(u) • n^- in (7.5) is not uniquely defined 
when replacing u by its approximation u p h £ V p {Vh)- Since u p h is discontinuous 
across element interfaces, q(u^) • n/< has two values on dK, one associated with 
u p h on the interior of element K and one associated with u p h on the exterior of the 
element. The ambiguity in evaluating the boundary flux in (7.5) is eliminated 
by replacing it with a numerical flux function: 


, / p,int(K) p,ext(K) 

n K,e\ U h > U h 


where u p ’ int{K) = u p h (x,y, •) ledge e of element k 

ext , K \ u h{ x iV-i Oledge e of element av fl T_ = 0 

u p h ' ex ' ' = < 

l d{t), di< n r_ / o 

The numerical flux function is required to have certain properties to 


enhance the convergence properties of the approximate solution: 
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1. h K t (u^’ tnt ^ K \ u£’ ex< ( A )) is a monotone flux function, i.e. h Kt is non- 

t . • p,int(K) i • p,ext{K) 

decreasing m u h v and non-increasing in ^ 

2. ^ Ke (v) is consistent, i.e., h Ke (u^u) = q( u ) • n A, e 

3. h K e (u% xnt ^ K \ v% ext ( K ^) is Lipschitz continuous 

4. h K < ,(u p h ’ tntf ' K \ u v ^ ext ^) is directionally consistent in the sense that 

, , p,int(K) p,ext(K)\ _ . , p,int(K c ) p,ext(K e K 

ri K,e\ U h i u h ) ”'I< e ,e\ U h 1 U h ) 

Note that this last property enforces continuity of the normal flux at element 
interfaces. Any function which satisfies the properties listed above can be used 
in this numerical scheme. Some examples of fluxes satisfying these conditions 
are the numerical fluxes of Godunov [23], Enquist-Osher [17], and Roe [46]. 

The semi-discrete problem which results from the discontinuous Galerkin 
approximation to (7.1) can now be written for a typical element K € Vh : 

Find u p h (x,t)\ K E Q p k[K) x C 1 ( 0,T] such that 


/ <(x,0)| /f u(x)dx = J u 0 (x)v(x)dx and 


d_ 

dt 


J uj(x, t)v(x) ix iy = / q(aj(x, ()) ' Vv(x)<ix 

eedK\an Je 

f h KA u h ,nt{h \ d ( t )) v ds 

e£d Kr\T~ Je 

/ q«) • n K.. u ds ( 7 - 6 ) 


eedK nr+“ 


for every v(x) € Q p x(K ) 
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Note that the integral over the element boundary is represented as the sum 
of integrals along element edges where n Ke denotes the outward unit normal 
to edge e of element I\ . For the linear conservation law, (7.6) reduces to the 
methods (2.54) and (2.55) provided that the numerical flux reduces to 



V,int(K) 

u h 


p,ext(K) 

i u h. 


+ <’ i,(K) )I3 ■ n 1(i , 




p,ext{K) p,int(I\) 


Wi 


•) 


when q(u) = /3u. 


7.2 Runge-Kutta Time Discretization 

The system of ordinary differential equations in (7.6) can be discretized in time 
using one of many time marching schemes, e.g. forward Euler, backward Euler, 
or some predictor-corrector method. These schemes combined with acceleration 
techniques, e.g. local time-stepping or multigrid, are good choices for problems 
where only the steady-state solution is of interest. For truly time-accurate 
solutions, however, it is necessary to have the order of accuracy of the time 
discretization to be equal to the order of accuracy of the spatial approximation. 
Using classical interpolation-type error estimates, Cockburn et. al. [12] show 
that the Galerkin spatial approximation with uniform p elements satisfies the 
estimate 

IIMO + V ' q«)lk~(Q) < C*<- +1 >|q(u)|„,„ +J ,„ (!1) (7.7) 

where we have used the abstract notation of Cockburn, Hou, and Shu [12] to 
represent (7.6), 


(7,8) 
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which formally requires the inversion of the ( p K + 1) x (p K + 1) mass matrix 
in (7.6). Runge-Kutta methods can be used to obtain a time marching scheme 
for which the truncation error is of order A <( p K +1 ) for u p h \ h . € Q P (K). 

Let At denote the time step increment and assume the approximate 
solution u p h (x, t n ), t n = nAt is known. We use a special class of explicit Runge- 
Kutta schemes written in the classical form as 

S 

u p h (x,t n+1 ) = u p h (x,t n ) + ^2b,k t 

»-l 

where k t = A tLh(u*(x, t n ) + ^ a i3 kj) (7.9) 

3 = i 

which is represented by the Butcher-array [9] 


0 




C 2 

a 21 



C s 

®sl 

^SfS — 1 



b\ 

- 

b s 


where C{ = £j=i a ij define the abscissa t{ — t n + C{At for a non-autonomous 
system which would result, for example, if the conservation law contained a 
time-dependent source term. The special class of schemes can be written in 
multi-stage form as 

m(0) = 

u (i) = n h ^[a tk uW + At/3 lk L h (uW)},i = \,...,s (7.10) 

k~0 

K(->in+ 1 ) = 

where and /?,*. are free parameters and 11^ denotes a projection operator 
to be discussed later. For consistency, > 0 and ^JZ.q = 1. We use the 
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1 1 
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3 6 


Table 7.1: Parameters for the TVBM Runge-Kutta scheme 


parameters listed in Table 7.1 which were derived by Shu [48] to result in a time 
discretization that is TVD (Total Variation Diminishing) when combined with 
finite difference schemes in one space dimension. For multiple space dimensions, 
the parameters listed in Table 7.1 result in a scheme which is TVBM (Total 
Variation Bounded in the Means) as defined in the following lemma. 

Lemma 9 (Cockburn, Hou, and Shu [12]) 

Let ao = inf un(x) 
xefi 
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b 0 = sup u 0 (x) 
xefi 

UK = f j^j j K U h{*iQ\K d x 

w{ d ^ f ul + 6 AtL h (u p h ) (7.11) 

where 8 > 0 is a parameter. Suppose that the operator Lh satisfies the following 
maximum principle: 

uk € [ a,b } — > wk € [a — Bh°) , 6 + £?/*"] (7-12) 

where a and B are non-negative parameters. Then the Runge-Kutta time dis- 
cretization defined by (7.10) is Total Variation Bounded in the Mean, i.e., 

Uf( E [«o — snBh 0 ) ,bo + snBh°) ] (7-13) 

Proof: By induction. See Cockburn, Hou, and Shu [12] for the case when 
a = 2. The case when a is arbitrary follows directly. H 

Remarks: 

(i) The second- and third-order schemes in Table 7.1 require little stor- 
age since at each stage there is only one non-zero parameter ft. 

(ii) The fourth-order scheme is the classical fourth-order Runge-Kutta 
method. For this and higher-order schemes, nearly all intermediate 
solutions must be stored since most of the parameters are non- zero. 


I 
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7.3 Enforcing the Maximum Principle 

A crucial step in the proof of the TVBM property is the assumption that the 
spatial operator satisfies the maximum principle (7.12). The local projection, 
U h in (7.10), is designed to enforce the maximum principle by modifying the 
solution at each stage in the Runge-Kutta scheme. Cockburn, Hou, and Shu [12] 
derive conditions on the solution on the boundary of an element so that (7.12) 
is satisfied. These condition are expressed as limits on the deviation of the 
solution on an element boundary from its mean values based on the difference 
between the element mean and the mean values in neighboring elements. The 
’’stencil” of elements used to define these limits includes neighboring elements 
as well as neighbors of neighboring elements. (See [12] or [6].) Note that (7.12) 
enforces the TVB property on solution mean values; it says nothing about 
the actual solution in an element. Moreover, values on the boundary of an 
element are insufficient to define a p-unisolvent element. Here we use a simple 
projection strategy which overcomes some of these difficulties. 

The overall strategy for the local projection is as follows: (1) identify 
elements in the neighborhood of a shock as those where the jump in the solution 
along the element boundary [u£] > (2) If p k > perform an intermediate 

projection P : Q P (K) — > Q l (K) (for p K = 1, P is the identity operator); (3) 
perform a local projection, II , on the result of step (2) to enforce the following 
condition: 

U K P{u p h \ K ) e [min(«A-; € dK), ma,x(u K ] u Ke , e € dK)] , Vx 6 I< (7.14) 

These conditions are the multi-dimensional extension of those used by Van 
Leer [51] in one space dimension and are identical to those used by Barth [5] 



in reconstructing piecewise linear polynomials from mean values in his higher- 
order finite volume schemes. 
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The limits in (7.14) can also be written in terms of the deviation of the 
solution from its mean value. Denoting the deviation from the mean by u k 
where u v h | = Uk + Uk, (7.14) becomes 

U f<P(u p h \ ) e [mm(0,u Kc - u K ),max(0,u Ke - S K )] , Vx € K (7.15) 

e(E.oh eEaft 

In order to maintain conservation, the projection must also preserve the mean 
value, 

Ri< p ( u h\ I( ) = «k (7.16) 

7.3.1 The Projection P : Q P {K) -> Q\K) 

There are several choices for the projection of a higher-order solution onto 
the space of bilinear shape functions; the only requirement is that it preserve 
the mean value. A computationally efficient approach for the projection P 
is to simply truncate the higher-order terms in polynomial representation of 
u v h \ k- To enforce the conservativity condition (7.16), the bilinear part of <\k is 
augmented by the mean value of the truncated higher-order terms. To describe 
this procedure, the approximate solution is decomposed into its linear and 
higher-order contributions: 

U h \ K — U linear T U h.o . 

4 (P K + 1) 2 

= E «.•&(*, y)+ H ai(j>i{x,y) (7.17) 

i=l t=5 


where <j>i(x,y ), i = 1,...,4 are the standard bilinear shape functions and 
u,-, i — 1,...,4 are the solution values at the nodes of the element. The 
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element mean value is thus 


V-K — ^ linear "4" ^h.o. 

4 _ ( P/f +l) 2 

y ' “h y ) ciifii 


(7.18) 


i= 1 i=5 

Let u? ew denote the degrees-of- freedom of P(uh\x), then its mean value 


is 


P«\ K ) = (7.19) 

t=l 

The conservativity condition (7.16) and the property of the bilinear shape func- 
tions that 4>% = 1 imply that 


i=l 


i= i 

4 4 

= + (H <f>i)uh. 0 . 

i= 1 «=1 

4 

^ ' (^i "t~ ^)i.o.)^i 
«=1 

= u t + u h . 0 .,i = 


(7.20) 


Combining the definition of u/i, 0 . in (7.18) with (7.20) results in the final 
definition for P(u p h | ): 


4 (P/c+l) 2 

p (“'L) - £(“( + S a i<h)M x <y) 


(7.21) 


i=l 


j-5 


7.3.2 The Local Projection IT/<- : Q*(K) — ► Q X (K) 

Construction of a projection which results in a function satisfying (7.15) for 
every point in the element is straightforward for linear elements since the max- 
imum solution values within the element occur at the nodes (vertices) of the 
element. 
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Let Pj, i = 1, • • • ,4 denote the degrees of freedom associated with the 
projected solution, 

u kP(uI\ k ) = Y,M x iV)Pi (7.22) 

i = 1 

which can be expressed in terms of the mean value and the deviation from the 


mean, 


n*p KL) = n K p(u;n + u K p(uiu 


(7.23) 


Evaluating (7.23) at the nodes yields a direct expression for the degrees 
of freedom of the projected solution, 


P , = n K p«| K ) + n, f PK! K )( I r,jf); > = i, ■■•,4 


(7.24) 


where (xf ,yh) ,i. = 1 , - • • , 4 are the coordinates of the element vertices. All 
that remains to define the local projection, 11^-, is to select values for 


n K PKi K )(*r,»r)< = !,•••, 4 

so that the monotonicity conditions (7.15) and the conservativity condition 
(7.16) are satisfied. To simplify the notation, let 


ft = ni<P(<\ K M,V?) = Tl K P(ul\ K )(x«,y!<)-n K P(ul\ K ) 

= «*(*?,»*) = «SL(*! t ’, »*)-«*■ , 1 = 1, •••.4 


u 


denote the deviation of the projected solution and the actual solution from the 
mean value, respectively. A strategy similar to Barth [5] is used to select p, so 
that (7.15) is satisfied: 


~ = f min(«„ : 
1 \ max(n,, 


min(u,-, max(«x; ux e , e G OK)) if Hi > 0 
min(u^; uk c , e G OK)) if u,- < 0 


(7.25) 


Note that (7.25) represents the minimum modification to u p h \x that satisfies 
(7.15) without regard for the conservativity condition. 
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The conservativity condition is enforced by appropriate scaling of pi. 
From (7.23) we have 


M<l t ) = 

1 = 1 

= J2M^K p ( u h\ K ) T Pi) from (7.24) (7.26) 

3 = 1 

and since IIa-P(«£| a . ) is a constant and 4>i — 1, (7.26) becomes 

£&W = ° (7.27) 

t=l 

To enforce conservativity, the appropriate p; are scaled by a factor a < 1 which 
is is a ratio of the positive and negative parts of (7.27). This scaling process is 
described by the following algorithm: 


Let S + = {i : pi > 0} 
S~ = {i : Pi < 0} 

p* = E ii ft 

t€5+ 

P- = E ft 

i&S~ 


If 

p- 

P < P + , then <j = •— and p; t-> erp, ; Vi € 5 + 


If 

p+ 

P > P + , then <7 = -p- and pi i— > crp, ; Vi G P 

(7.28) 

If 

P~ = 0 or P + = 0, then p, = 0, i = 1, . . . ,4 



This completes the definition of the projection for p — 1 elements since from 
(7.24), (7.16), and (7.23) we get 

n*p«|J = + p«' 


(7.29) 
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7.4 Numerical Results for Burger’s Equation 

The performance of the discontinuous Galerkin method is investigated by solv- 
ing the two-dimensional Burger’s equation with periodic initial and boundary 
conditions: 


ut + f(u)x + g(u) y 
where f(u) 
u(x,y,0) 

«( i,y) 
«(*»!) 


o (x,y) € (-1,1) X (-1,1), t > 0 

u 2 

»W = T 

11 ., 

- + -smTT {x + y) 

«(-i ,y) 

u(x , —1) 


(7.30) 


The solution to (7.30) is smooth until the time t = 0.45 when two shocks 
form diagonally to the domain boundaries. The exact solution at t = 0.1 and 
t = 0.45 are shown in Figure 7.1 where a 40 x 40 mesh is used to generate 
the contours and 100 points are used to generate the distributions along the 
domain diagonal. The exact solution for t > 0 is computed using the method 
of characteristics. 


The time accurate solutions presented in this chapter were obtained 
using the Godunov numerical flux function [23] 


f min„ £ [ UltU2 ] q(u ) • n X e for «i < u 2 

\ max ue [„ 2iUl ] q(u) ■ n Ke for Ui > u 2 


(7.31) 


and the Runge-Kutta time marching schemes defined in Table 7.1 using p + 1 
intermediate steps for elements with a polynomial approximation of degree p. 
To assess the accuracy of the method, the error in the approximate solution, 
e = u — u p h , is computed in the L 1 and the L°° norms. The L 1 norm is a global 
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measure of the error and is defined as 

||u — UaIUiCO) = \u-u p h \dxdy 

J Q 

= I \ u ~ u h\ dx d v (7.32) 

Kev h JK 

The element integrals in (7.32) are evaluated using numerical quadrature with 
nine integration points in each local coordinate direction to accurately compute 
the error near discontinuities. The L°° norm is a local measure of the error and 
is defined as 

ll«-«fcllL~(0) = max . |u-uj| 

— max 
K&> h 

w max (7.33) 

{xi,yj)eK 

where (a:,-, t/j), i,j = are the coordinates of the integration points in 

element K . 

Numerical results are shown in the form of distributions of the solution 
along the domain diagonal. The actual discontinuous solution is displayed by 
subdividing each element into 10 smaller elements, evaluating the approximate 
solution at the vertices of the sub-elements, and assuming a linear distribution 
of the solution in each sub-element. Discontinuities in the solution at element 
interfaces appear as vertical lines connecting two circles in the distribution 
plots. These circles represent the solution values at the vertices of an element. 
Since the solution is periodic, only one period is shown in the distribution plots. 

To assess the accuracy of the Runge-Kutta discontinuous Galerkin method 
without the local projection, results are presented at t = 0.1 when the solu- 
tion to (7.30) is smooth. The solution obtained on an 8 x 8 mesh of uniform 


max It/ — u v Ak 

{x, y )eK 
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p — elements (p = 1, • • • , 4) is shown in Figure 7.2. Note that the approximate 
solution is discontinuous at element interfaces even though the exact solution 
is smooth. These discontinuites are small, 0(/d p+1 )), and thus decrease as the 
mesh is refined or as p K is increased. 

The rate of convergence of the error in the Zd-norm for a sequence of 
uniform mesh refinements is shown graphically in Figure 7.3. Each line on the 
graph represents a sequence of meshes with a fixed element polynomial degree 
p. The slope of the lines on the log-log scale is nearly p + 1 indicating that the 
error decreases at the rate of h p+1 as the mesh size h decreases. In other words, 
the method is (p+ l)-order accurate. 

A graph of the error as a function of the total number of unknowns 
in the solution, Figure 7.4, shows that the number of unknowns required to 
achieve a certain level of accuracy decreases by an order of magnitude for each 
increase in p. Figure 7.4 also shows that for a fixed problem size, more accurate 
solutions are obtained using higher-order elements. 

As the solution evolves, the local projection is required to control oscil- 
lations in the solution. To demonstrate the effectiveness of the projection at 
controlling oscillations, the solution obtained at f = 0.45 on a uniform 8x8 
and 40 x 40 element mesh with p — 1 and p — 2 is shown in Figure 7.5. The 
solutions obtained for p = 3 elements in indistinguishable from the p = 2 so- 
lution, and is therefore, not shown. In all cases, the projection is extremely 
effective at controlling oscillations at the shock. 

The error in the solution in smooth regions, computed on a subdomain 
which excludes the shock, is summarized in Figure 7.6. Each line on the graph 
represents a sequence of meshes with uniform p-elements. Note that the slope 
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Figure 7.6: Rate of convergence of the error in smooth regions at t = 0.45. 

(ft, = [-0.2, 0.4] x [-0.2, 0.4]) 

of the lines is p + 1 indicating that the order of accuracy of the method is 
preserved in smooth regions. 



Chapter 8 


Concluding Remarks 


8.1 Summary 

The development of a parallel ftp-adaptive discontinuous Galerkin method for 
hyperbolic conservation laws is presented in this work. A priori error estimates 
are derived for a model class of linear hyperbolic conservation laws. These 
estimates are obtained using a new mesh-dependent norm that reflects the 
dependence of the approximate solution on the local element size and the local 
order of approximation. The results generalize and extend previous results 
on mesh-dependent norms to ftp-version finite elements and to discontinuous 
Galerkin methods. 

A posteriori error estimates which provide bounds on the actual error 
are developed in this work. The a priori and a posteriori estimates play an 
essential role in the development of an ftp- adaptive strategy designed to deliver 
solutions to a specified error level in an efficient way. A generalization of 
the three-step ftp-adaptive strategy is developed using the error estimates to 
provide for detection of discontinuities in the solution and local ft-adaptivity 
when appropriate. 

Numerical experiments verify the a priori estimates and demonstrate the 
effectiveness of the a posteriori estimates in providing reliable estimates of the 
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actual error in the numerical solution. The numerical examples also illustrate 
the ability of the /^-adaptive strategy to provide super-linear convergence rates 
with respect to the number of unknowns in the problem, even for discontinuous 
solutions. 

A parallel implementation of the discontinuous Galerkin methods is pre- 
sented which takes full advantage of the local character of the method and 
results in nearly optimal speedups on hp- meshes. 

Extensions of the discontinuous Galerkin methods to nonlinear hyper 
bolic conservation laws are also presented. Numerical results illustrate the 
effectiveness of the method at delivering high-order accuracy in smooth re- 
gions. A local projection, designed to control nonlinear stability, is shown to 
eliminate oscillations and provide high resolution of discontinuities. 

8.2 Conclusions and Future Work 

The study reported in this dissertation represents a significant departure from 
conventional studies on the numerical solution of hyperbolic problems. The aim 
was to produce new schemes which deliver very high accuracies, were readily 
parallelizable, were based on rigorous mathematical foundations, and which 
were capable of delivering exponential rates of convergence. All of these goals 
were accomplished for a model class of linear hyperbolic conservation laws in 
two-dimensions, and encouraging results were obtained on extensions to model 
nonlinear problems. 

Among the major conclusions drawn from this work are the following: 


1. The machinery of elliptic approximation theory can be extended to hp- 



122 


finite element approximations of hyperbolic equations using the notion of 
discontinuous Galerkin methods; this is made possible by the introduction 
of special bilinear and linear forms which depend upon mesh parameters, 
the mesh size h K of a cell K in the mesh and the spectral order p of the 
shape functions characterizing local approximations over the cell. 

2. The use of the new mesh-dependent norms makes it possible to derive, for 
the first time, a priori error estimates for non-uniform hp - approximations 
of linear hyperbolic problems; these estimates are quasi-optimal, and the 
estimated rates of convergence are fully confirmed by numerous numerical 
experiments. 

3. Exponential and/or super-linear rates of convergence are obtained, even 
for solutions with very low regularity; this justifies the use of /ip-methods 
and demonstrates their superiority over conventional methods for a model 
class of problems. 

4. Rigorous a posteriori error estimates are derived using a new version of 
the element residual method; these estimates provide computable mea- 
sures of local (elementwise) error with remarkable accuracy and provide 
a reasonable basis for assessing solution quality and for adaptivity. 

5. Equipped with both a priori and a posteriori estimates; an effective hp- 
adaptive strategy is developed which can be parallelized and generally 
gives a good hp - mesh in three or four steps; this work represents an 
extension and generalization of the three-step scheme for non-uniform 
hp- meshes; it exploits the unique feature of the a priori estimates for 



hyperbolic problems, in particular, the loss of the rate of convergence in 
the vicinity of a discontinuity. 
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6. The adaptive scheme developed here is quite robust and effective for the 
example problems tested; results suggest that it is possible to specify 
target global error in an appropriate norm and to then reach that error 
quite accurately in three or four steps. 

7. New versions of load-balancing schemes based on recursive bisection pro- 
vide for a domain decomposition well-suited for /ip-version disontinuous 
Galerkin methods; the schemes, which are in an early stage of develop- 
ment, still provide a reasonably balanced computation when implemented 
on a 16-processor Intel iPSC 860 computer. Near linear speedups were 
realized on a model problem. 

8. Extensions of the methodologies to nonlinear problems appear to be pos- 
sible; preliminary studies suggest that with a proper projection to control 
oscillations near discontinuities, very high accuracies can be obtained with 
hp - schemes using the discontinuous Galerkin method; while much work 
remains to be done in this area, the high convergence rates and accura- 
cies observed in the numerical experiments on nonlinear problems suggest 
that further studies in this area may be very fruitful. 
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